Il me semble que la thèse s'apparente à un voyage, du moins tel que le décrit Nicolas Bouvier dans 'L'usage du monde' : "Un voyage se passe de motifs. Il ne tarde pas à prouver qu'il se suffit à lui-même. On croit qu'on va faire un voyage, mais bientôt c'est le voyage qui vous fait ou vous défait". Aux antipodes de l'aventure solitaire, cette thèse s'est révélée riche en rencontres. Des rencontres auxquelles je souhaite rendre ici toute leur importance. Et puisque la mode est à l'information compressée ;-), tâchons de faire court et bien ! C'est à mon encadrement pendant ces trois années que vont naturellement mes premières pensées. Merci à toi Cyril, l'instigateur de cette thèse, pour avoir partagé avec moi ton expertise et ta vision de la reconstruction tomographique. Merci également pour ton soutien sans faille et tes qualités humaines. C'était un réel plaisir de travailler à tes côtés ! Merci à toi Nikos, pour avoir toujours soutenu mon ambition en m'incitant à viser les étoiles. Elles ne sont pas encore atteintes mais restent dans l'objectif. Merci à toi Arthur, pour la pertinence de tes apports dans un domaine qui ne t'était pas familier, toi le 'statisticien'. Merci aussi pour ta bonne humeur et la narration d'une certaine anecdote -que je ne relaterai pas ici, non-mais qui a beaucoup amusé l'auditoire de ma soutenance... Merci à toi Yves, pour avoir contribué à développer ma curiosité clinique et m'avoir appris que j'écrivais comme Kirkegaarde -ou peu s'en faut. Merci aussi pour ces innombrables anecdotes culturelles dont tu as le secret et qui rendent ta conversation toujours extrêmement stimulante. Merci à vous Gilles et Elisabeth, pour votre regard critique, votre soutien dans les grandes échéances et votre bonne humeur au quotidien.

Ma reconnaissance va également à ces experts reconnus, du monde de la reconstruction, de l'imagerie ou encore de l'optimisation, et qui m'ont fait l'honneur d'évaluer ma recherche à l'aune de la communauté scientifique. Merci à vous Bruno, Alfred, Jean Christophe et Jerry, pour votre sérieux et l'intérêt que vous avez porté à mes travaux. Votre lecture critique et les discussions qui en ont découlé ont permis de mettre en exergue les faiblesses mais aussi les potentialités de cette recherche, un regard neuf et objectif pour lui donner encore plus de relief.

Le quotidien aurait été parfois bien morne sans les collègues des différentes entités avec lesquelles j'ai interagi. Un grand merci donc à vous collègues de GE Healthcare, pour m'avoir intégrée à votre équipe et associée à vos activités de recherche. Merci aussi pour tous ces bons moments en votre compagnie que ce soit permanent ou stagiaire, à l'occasion d'une pause café, d'un repas, d'un trajet (en bus, voiture, ou avion), d'un footing et j'en oublie... Grâce à vous, ma vie est désormaix plus caféinée (ou devrais-je dire, cappuccinée ?). Merci à vous collègues de Supélec, toujours aussi accueillants malgré mes allers et venues en 'courant d'air'. Merci plus particulièrement à toi Luc, pour ton temps, ta générosité et toutes ces grandes discussions sportives. Merci à vous collègues de courte durée de UPenn. Merci Christos pour m'avoir si gentiment intégrée à l'équipe. Vous avez contribué à me faire dépasser les clichés, à me faire aimer Philadelphie et cet 'autre côté' de l'Atlantique. Pour peut être un jour y retourner ? Mes derniers remerciements s'éloignent du cadre purement professionnel pour toucher un cercle plus personnel. Merci à toi Vincent, que je retrouve maintenant avec plaisir comme collègue. Sans ton soutien je n'aurais sans doute jamais osé me lancer dans cette aventure qu'est la thèse. Merci à vous Thomas et Aymeric, mes deux folles, pour votre soutien indéfectible, votre bonne humeur communicative et votre incroyable capacité à me faire rire, même face à l'adversité. Aussi parce que trouver une pile 9V à 30 minutes du début de ma soutenance, ça n'a pas de prix ! Merci à vous Charlotte et Dustin, pour compléter naturellement ce quatuor de choc qui m'a donné tellement d'énergie durant mes derniers mois à GE. Merci à vous amis rencontrés au détour des conférences. Nous nous reverrons sans doute un jour, à l'autre bout du monde. Merci à vous athlètes et entraîneur APY-Bucois, pour tous ces moments partagés sur le stade. Vous avez été et continuez d'être une délicieuse bouffée d'oxygène. Merci à vous amis de toujours, pour avoir accepté mes choix et les avoir toujours soutenus, malgré mon peu de disponibilité et le stress engendré. Merci à vous papa et maman, pour cet amour de la connaissance que vous m'avez transmis et ces ailes que vous m'avez données. Enfin, merci à toi, Nadège, petite soeur chérie, pour être telle que tu es. Continue de me surprendre. Je te souhaite d'également beaux voyages.

This work addresses the reconstruction of X-ray tomographic acquisitions in the clinical context of interventional procedures. Interventional radiology refers to minimally-invasive procedures, that consist in inserting and manipulating medical tools under image guidance through the vascular system of the patient to the pathology location, where the tools are used to deliver the treatment. Image guidance of interventional procedures was historically carried out through two-dimensional real-time projective images acquired in low-dose X-rays with C-arm systems. Current C-arm systems present both the advantages of a tomographic tool -since the rotation of the C-arm gantry around the patient allows for collecting a data set of X-ray measurements for a large amount of angles-and of a facilitated patient access. Because the clinical routine commonly involves enhancement of the blood vessels through intraarterial injection of a contrast agent, rotational angiography provides a three-dimensional visualization of the vasculature, that has proven to be of great help for improving the management of the medical gesture.

Tomographic reconstruction is computed analytically through filtered backprojection whose discretization defines sampling requirements that are central to this work. In particular, a rapid acquisition is needed to guarantee data consistency between the projection data and hence, feasibility of the threedimensional reconstruction. A sampling trade-off is defined for rotational angiography, that is not as favourable as the one defined for diagnostic imaging with CT scanners because of the low framerate of the detector (responsible for uniform subsampling) and the mechanical constraints that are inherent to the current C-arm systems (rotation speed not high enough to avoid patient motion, which results in limited range subsampling), thus remaining the main limitations to more use of tomography in image-guided therapy.

Although the visualization of the densest structures (e.g. bones, contrast-enhanced vessels) is sufficient for the current clinical practice with the framerate delivered by C-arm systems, efforts for the visualization of the soft tissues are underway. In addition, the angiographic data are frequently corrupted by physiological motions (e.g. cardiac motion, respiratory motion) resulting in reconstruction artifacts that may hamper the medical interpretation. Considering more specifically the case of cyclic motions, a widespread approach in the diagnostic imaging of dynamic data with CT scanners consists in performing multiple rotations and selecting one subset of consistent views corresponding to the same phase of the motion cycle, for which filtered backprojection reconstruction is feasible. With C-arm systems, however, one single rotation is performed only, resulting in poorly sampled phases.

In this work, we rely on compressed sensing [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF] -a recent mathematical theory that makes more efficient use of the collected samples and that met considerable impact in the statistics and signal communities-to propose original iterative reconstructions that are adapted to rotational angiography. The compressed sensing theory defines assumptions under which the recovery of subsampled data is possible through the minimization of a least-square data fidelity term that is combined with sparse constraints expressed as ℓ 1 -penalties. Although the compressed sensing assumptions are not strictly valid for X-ray tomography, the theory was illustrated with exact reconstruction of a piecewise constant phantom from poorly angularly sampled data that simulated the X-ray sampling process. This opened tremendous perspectives for X-ray tomography. Numerous reconstruction algorithms based on compressed sensing have been recently proposed in the literature. To this purpose, three-dimensional reconstruction is cast as an optimization problem that is handled by transform-domain sparsity wherein sparser signals allow CHAPTER 1. INTRODUCTION for significant mitigation of the subsampling artifacts compared to conventional analytical reconstruction techniques. Thus the compressed sensing reconstruction of the coronary arteries based on the minimization of the image ℓ 1 -norm was proposed for cardiac angiography [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF]. Motivated by the desire to reduce dose, iterative reconstruction based on total variation minimization was proposed for tomosynthesis [START_REF] Sidky | Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved imagereconstruction algorithms[END_REF]. To reconstruct dynamic data, a particularly attractive approach is the PICCS algorithm [START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF] that uses a sparse constraint based on a prior image -that may be the static filtered backprojection reconstruction-, so that the reconstruction, schematically, determines motionless areas from the complete acquisition, while motion-blurred areas are determined from reconstruction of well-identified subsampled motion phases only. This strategy, however, is not applicable to intraarterial contrast injection.

In the present work, we propose an application of the compressed sensing theory and a generalization of the PICCS algorithm that are adapted to rotational angiography with intraarterial injection. An iterative implementation based on filtered backprojection (iFBP) is developed, that allows for efficient computation of the least-square solution. The algorithmic framework relies on an iterative proximal splitting scheme [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF], in which the minimization of the least-square fidelity term of iFBP is alternated with a proximal operator that allows for incorporating a series of ℓ 1 or ℓ 2 penalties into the reconstruction process. This approach is derived for single and multiple sparse penalties. Observing that compressed sensing is not reduced to ℓ 1 -norm minimization, a homotopy continuation strategy is investigated, that presents relationship with orthogonal matching pursuit [START_REF] Donoho | Sparse Solution of Underdetermined Systems of Linear Equations by Stagewise Orthogonal Matching Pursuit[END_REF].

The evaluation of the compressed sensing reconstructions developed in this work is performed on both simulated and clinical data with focus on injected datasets: severe streak artifacts may arise, that originate from subsampled temporal variations of the contrast-enhanced vessels. The application to a selection of clinical abdominal and cerebral datasets suggests a domain of applicability:

• clinical cases of subtracted reconstruction, contrast-enhanced reconstruction, and reconstruction with a limited range, for which a static reconstruction is proposed;

• clinical cases of digital subtraction rotational angiography, small pulsatile displacements of the contrast-enhanced vessels and contrast variations, for which a dynamic reconstruction is proposed.

The reconstruction results outperform current analytical reconstruction and clearly demonstrate the interest of using compressed sensing with homotopy in rotational angiography.

Publications

Most relevant results of this work can be found in the following publications:

Conference publications (with review committee)

• H. Langet Chapter 2

X-Ray Computed Tomography

X-ray computed tomography (CT) is a technique that measures the density distribution of the human body through X-ray measurements. To this aim, the body is virtually carved into thin two-dimensional (2D) cross-sectional slices1 . The rotation of the imaging system around the patient allows for collecting a data set of X-ray measurements for a large amount of angles, from which recovery of the density distribution is feasible. Tomographic reconstruction is an inverse problem for which an analytical formulation exists through the inversion of the 2D Radon transform, that is presented in Sec. 2.1 together with its discretization and associated sampling requirements that are central to this work. Since early experimental X-ray imaging (Roentgen, 1895), radiological systems have benefited from steady technological progresses. We introduce in Sec. 2.2 the physics underlying X-ray interactions with biological tissues. X-ray detection for diagnostic imaging is presented in Sec. 2.3 through the technical characteristics of the X-ray CT scanner.

Mathematical basis of tomography

The purpose of this section is to introduce the mathematical concepts and tools that underlie tomographic reconstruction. Let us describe an image as the function of two spatial variables describing some physical property of an object. Computed tomography is an indirect imaging technique: there is access to the function through its integrals over straight lines. An inverse operator is therefore required to determine the image from the set of integrals. This is achieved by the inversion of the Radon transform in parallel-beam geometry. For more details on the concepts and mathematical demonstrations we refer the interested readers to [START_REF] Herman | Image Reconstruction from Projections. The Fundamentals of Computerized Tomography[END_REF][START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF].

Radon transform in two dimensions

Let us consider a 2D Cartesian coordinate system where O denotes the origin, -→ x the positive x-axis unit vector and -→ y the positive y-axis unit vector. Let us denote f a function where f (M ) = f (x, y)

is the value of some physical property at point M of coordinates (x, y). It has a finite support Ω, i.e. f (x, y) = 0 for points (x, y) that are located outside the reconstruction field of view, that is usually a circle of diameter D.

We introduce the coordinate system ( -→ u , -→ v ) that is defined by rotating ( -→ x , -→ y ) at angle φ with respect to the positive x-axis, such that u = x cos φ + y sin φ v = -x sin φ + y cos φ .

As illustrated in Fig. 2.1, we define a parallel-beam geometry as the set of lines D PB (u, φ) that are parametrized by:

u ∈ ] -∞, +∞[ φ ∈ [0, π[ .
The Radon transform of function f , that is denoted [Rf ], is defined as the complete collection of line integrals over the parallel-beam geometry [START_REF] Radon | Über die Bestimmung von Funktionen durch ihre integralwerte langs gewisser Mannigfaltigkeiten[END_REF]. The projection of f at angle φ, that is denoted p φ (u), is thus expressed as

p P B,φ (u) = [Rf ](u, φ) = M ∈DPB(u,φ) f (M )dM = +∞ -∞
f (u cos φ -v sin φ, u sin φ + v cos φ)dv

(2.1)

The projection of point (x 0 , y 0 ) is described by the parametric curve u(φ) = x 0 cos φ + y 0 sin φ, i.e. a sinusoidal curve in the (u, φ)-diagram. This diagram is thus called sinogram, that is the standard way of displaying the Radon transform of an image. Figure 2.2 shows an image and its sinogram: the projection data consist of the sampled superimposition of all sinusoids, each one weighted by value f (x, y).

Radon transform inversion

The inversion of the Radon transform is easily understood in the context of Fourier analysis.

Projection-slice theorem

Let us define F 1D Rf the one-dimensional (1D) Fourier transform of the Radon transform along u as:

[F 1D Rf ](k u , φ) = +∞ -∞ p P B,φ (u)e -i2πkuu du
Let us now introduce F 2D f the 2D Fourier transform of the function f by writing:

[F 2D f ](k x , k y ) = +∞ -∞ +∞ -∞
f (x, y)e -i2π(kxx+kyy) dxdy The projection-slice theorem states that the 1D Fourier transform of [Rf ] is equal to the 2D Fourier transform of f according to

[F 1D Rf ](k u , φ) = [F 2D f ](k x , k y ) with k x = k u cos φ k y = k u sin φ (2.2)
It follows immediately from Eq. (2.2) that there is a one-to-one correspondence between the Radon transform and the 2D Fourier transform and that the Radon transform completely describes any function f through the polar sampling of its 2D Fourier transform.

Filtered backprojection (FBP)

Function f can be recovered from its Fourier transform by the inverse Fourier transform: Relying on the projection-slice theorem given in Eq. (2.2), [F 2D f ](k u cos φ, k u sin φ) can be replaced with [F 1D Rf ](k u , φ). Finally, by symmetry of trigonometric functions and symmetry of the Radon transform the integration in φ can be limited to half of the angular domain. Provided that u = x cos φ + y sin φ, we obtain:

f (x, y) = +∞ -∞ +∞ -∞ [F 2D f ](k x ,
f (x, y) =
where h(u) is the ramp filter defined by the frequency response |k u |. This inversion formula has thus two steps:

• a 1D filtering of the Radon transform with the so-called ramp filter;

• a simple integration over φ that is called backprojection and corresponds to the adjoint operator of the Radon transform.

Discretization and sampling requirements

Digital systems are only able to collect discrete measurements. Up to now, we assumed for the simplicity of the mathematical analysis that continuous measurements are gathered and we presented an integralform solution. The numerical implementation of FBP however requires the discretization of the inversion formula.

Shannon-Nyquist criterion

Mathematically, a 2D function that is spatially sampled at rate 1 ∆x along the x-axis and 1 ∆y along the y-axis is modelled as the product of the continuous function and the 2D Dirac comb δ [START_REF] Roubine | Distributions-signal, ser. Collection[END_REF] of spatial periods ∆x and ∆y:

f s (x, y) = f (x, y) +∞ i=-∞ +∞ j=-∞
δ(x -i∆x, y -j∆y).

As a result, the spectrum of f s (x, y) is a periodisation of the spectrum of f (x, y): it consists of the original spectrum of which have been added its replicas at sampling frequency and associated harmonics.

If the shifted spectral patterns overlap, in a process called aliasing, they become indistinguishable from one another. To recover an estimate of the measured signal that is as close as possible to the original one, none of the replicated frequencies should overlap. Let us assume f is band-limited of maximum frequencies k max x and k max y . All else being equal, the well-known Shannon-Nyquist sampling criterion states that, to avoid frequency aliasing, it is desirable to gather the measurements at a rate above the critical shift that is defined as twice the spectrum upper bound:

1 ∆x ≥ 2k max x 1 ∆y ≥ 2k max y . (2.4)
Data collection at a rate below the Shannon-Nyquist criterion, where aliasing occurs, is referred to as subsampling.

Sampling requirements

In tomography, the sampling of the Fourier domain is on a polar grid instead of a Cartesian grid. A polar coordinate system can be described both in terms of radial and tangential directions:

• The radial direction defines the number of measurements (or bins) that are collected per angle. The sampling requirements for the radial direction define the spatial sampling and can be assessed following the Shannon-Nyquist theory;

• The tangential direction defines the number of angles (or integral lines) that are collected. The sampling requirements for the tangential direction define the angular sampling. We shall see that the Shannon-Nyquist theory is not as helpful for determining the required number of angles.

Radial sampling requirement Function f has a finite support, and hence its spectrum has no upper bound. However, it is not the full spectrum that is measured, but instead its convolution with a rectangular window that models the data collection system. Thus, the maximum achievable spatial frequency k max u along u can be shown to depend on the resolution ∆u of the data collection system. Let us approximate the response of a measurement bin by the rectangular window function 1 ∆u rect ∆u :

1 ∆u rect ∆u (u) =    1 if u ∈ -∆u 2 , ∆u 2 0 otherwise
. The measurement at bin u, denoted p φ (u), corresponds to the convolution of the Radon transform p φ with the window function 1 ∆u rect ∆u , that is evaluated at position u∆u:

p φ (u) = 1 ∆u rect ∆u * p φ (v) v=u∆u
It follows that the magnitude of the spectral pattern that is replicated by sampling is attenuated by the sinc function:

[Fp φ ](k u ) = 1 ∆u k sinc π ξ - k ∆u ∆u • [Fp φ ] ξ - k ∆u ξ= k u ∆u
, and consequently, the maximum achievable frequency can be assumed to be the first zero-crossing point of the sinc function. Since the sampling in k u is polar, the sampling in k x (φ = 0 • ) and k y (φ = 90 • ) admits the same upper-bound. It follows from Eq. (2.4) that ∆x = ∆y = ∆u.

(2.5)

Hadamard introduced the concept of ill-posed problem to define problems that result in unstable solutions. A problem is said to be well-posed if 1. its solution exists, 2. its solution is unique, 3. it is well-conditioned, i.e. small differences in the projection data yield slight differences in the solution, while a problem that is not well-posed is said to be ill-posed. Tomographic reconstruction implies frequency amplification by the ramp filter. As ∆u tends to zero to sample higher frequencies and increase the spatial resolution of the discretized solution, the spectral amplification tends to +∞, thus violating Hadamard's property of a well-conditioned problem (point 3). In practice, a low-pass filter is used to band-limit the ramp filtering. The frequency representation of the ramp filter with a rectangular window function is provided in Fig. 2.4(a). The rectangular window provides relatively high resolution but also noise amplification, and introduces the so-called Gibbs oscillations. Both issues are classically handled with apodization functions (e.g. Hanning, Hamming, Shepp-Logan, Butterworth, • • • ) that provide a continuous transition between the ramp amplification and the band-limited windowing. Angular sampling requirement Let us assume that the tangential direction of the 2D Fourier domain is sampled by N integral lines that are uniformly spaced at constant angular interval ∆φ over 180 • . Considering the star pattern formed by all the integral lines that go through the origin of the image, we observe that the density of the tangential sampling decreases with distance from center, and hence is non-uniform. To avoid aliasing, it would be desirable to consider a maximal angular interval so that, in particular, the Shannon-Nyquist sampling criterion is satisfied for adjacent samples at the periphery of the support of f : ∆φ = 2∆u D .

It follows that

N = πD 2∆u . (2.6)
The angular patterns that are associated to full sampling and to uniform subsampling are schematically represented in Fig. 2.5(a) and in Fig. 2 sampling criterion exposed in Eq. (2.6) is not used as it may lead to an unnecessary high number of angles. To illustrate why it is so, Fig. 2.6 compares the reconstruction obtained from 64 lines and 256 presents many high resolution details despite the fact that it was obtained with an angular sampling far less favourable than the theoretical criterion, while a severe loss of spectral content is observed in Fig. 2.6(c), that cannot be compensated by the high angular sampling. The number of lines N is rather determined empirically and depends on the application: CT scanners thus use about 500 lines that are decreased to 150 lines with C-arm systems. Filtered backprojection yields an artifact-free reconstruction only if the set of measurements is infinite. To illustrate the effect of angular subsampling on the quality of the reconstruction, we simulated the acquisition of subsampled projection data from the imaged object that was displayed in Fig. 2.6(a). We sampled the image spectrum with 32 equiangularly spaced integral lines over 180 • . As a result, artifacts that typically consist of a set of intense straight lines become clearly visible with impact over the entire image, as shown in Fig. 2.7(a). These streaks emanate from high-contrast structures and go across the full image with angular directions that correspond to missing angles. The streak pattern is emphasized in Fig. 2.7(b) that corresponds to the difference between the Fig. 2.7(a) and Fig. 2.6(a). The impulse response of the band-limited ramp filter that is illustrated in Fig. 2.4(b) presents negative components around the central positive peaks. When the set of measurements is infinite, the positive and negative contributions compensate each other during backprojection. When the number of integral lines is decreased, non-compensation of the positive and negative contributions results in non-localized streak artifacts.

Angular range requirements

So far, we explored the sampling requirements assuming a rotation over 180 • . In this case, the discretized backprojection formula is:

f (x i , y i ) = π N n pPB,φn (x i cos φ n + y i sin φ n ). (2.7)
Let us now consider the case where the object is not static during the 180 • rotation so that the projections do not correspond to the same object depending on the angle and are thus called 'inconsistent'. If we sort the projections into M subsets E m of consistent projections, each subset defines a new angular

f m (x i , y i ) = π N n∈Em pPB,φn (x i cos φ n + y i sin φ n ), (2.8)
so that:

f (x i , y i ) = m f m (x i , y i ).
(2.9)

Note that we could as well define:

f m (x i , y i ) = π Card {E m } n∈Em pPB,φn (x i cos φ n + y i sin φ n ), (2.10)
where Card {E m } refers to the cardinal of E m , so that:

f (x i , y i ) = π N m Card {E m } f m (x i , y i ).
(2.11)

In the following, we quickly show that, with filtered backprojection (FBP), it is not advantageous to reconstruct consistent data over a limited range, but instead that less artifacts are seen when reconstructing inconsistent data over the full 180 • range. We simulated the acquisition in parallel-beam geometry of 300 equiangularly spaced integral lines over 180 • of an object that is seen in two different phases during the scan. The object is made of a moving structure of contrast 5000 over a static background of value 1000. Phase 1 that is shown in Fig. 2.9(a) is imaged from 0 to 120 • , while phase 2 that is shown in Fig. 2.9(b) is imaged from 120 to 180 • . The partial reconstructions for phase 1 that is shown in Fig. 2.9(c) and phase 2 that is shown Fig. 2.9(d) have all their structures severely affected by the restricted angular coverage. We observe empirically that it may be desirable to neglect the phase dependency, and thus implicitly assume that the object of interest remains stationary during the acquisition of the measurements and apply classical FBP. Since backprojection is an integration process, static FBP reconstruction that is given in Fig. 2.9(e) corresponds to the sum of all the partial reconstructions. The positive and negative contributions among neighbouring lines are combined, and no artifacts appear in the final reconstruction. When inconsistencies occur between the measurements, the reconstruction process is no longer able to properly combine the positive and negative contributions and artifacts result. Thus, the sampling artifacts of the static background in phase 1 are compensated by the complementary artifacts in phase 2, while the sampling artifacts of the moving structure do not compensate each other so that the entire image is deteriorated. In the following, we refer to these artifacts as motion artifacts. To study the dependency of motion artifacts with respect to the imaged object and the angular sampling pattern, we simulated the acquisition in parallel-beam geometry of 300 equiangularly spaced integral lines over 180 • of measurements for two two-phase objects with a static background of value 1000:

• the object that is displayed in Fig. 2.10(a) contains a large structure of contrast 100, that is affected by a uniform periodic displacement during data collection as illustrated in Fig. 2.10(b);

• the object that is displayed in Fig. 2.11(a) contains a small structure of contrast 5000, that is affected by either a uniform periodic displacement in Fig. 2.11(b), a nonuniform periodic displacement in Fig. 2.11(c) or a sudden contrast change in Fig. 2.11(d) during data collection.

When the magnitude of the temporal variation (i.e. intensity change) is low, the inconsistencies within the scan are small and yield a reconstruction that is deteriorated by localized mild blurring as shown in Fig. 2.10(c). When the magnitude of the temporal variation is large, the larger inconsistencies within the scan yield a reconstruction that is deteriorated by strong streaks as illustrated in Fig. 2.11(e), Fig. 2.11(f) and Fig. 2.11(g). With regard to Fig. 2.11(e), Fig. 2.11(f) shows that the severity of the spectral subsampling strongly depends on the angular sampling pattern associated to the time-points.

In particular, nonuniform subsampling results in more prominent artifacts than uniform subsampling.

Formation of the X-ray measurements

We have shown how to determine the 2D function of a physical property from its integral measurements. We now focus on attenuation of biological tissues to X-rays that is the property measured by the medical devices we shall consider further on. This section aims at providing some elements of the X-ray physics, that are necessary to understand the image formation process. We refer the interested readers to [START_REF] Dendy | Physics for Diagnostic Radiology[END_REF] for more details on the underlying physical principles of X-rays.

Production of X-rays

As schematically illustrated in Fig. 2.12(a), the main components of a Coolidge X-ray tube (Coolidge, 1913) comprise:

• an evacuated chamber that maintains a high vacuum and across which a high voltage (in kV) can be applied. X-rays are emitted out of the chamber through a thin window;

• a cathode that is housed in the evacuated chamber and that emits thermally-excited electrons by heating a metal filament (e.g. tungsten);

• a metal anode (also called target) that is housed in the evacuated chamber and that converts electrons into X-ray photons, while being able to quickly dissipate the large amount of generated heat. The material and the design of the anode are essential features for the tube efficiency and lifespan;

• an envelope that houses and protects all parts. The common materials used for manufacturing for envelopes are glass, metal or ceramic. A displacement of charges from the cathode to the metal anode is induced through the application of a difference of potential that can go up to 150 kV. The electrical charge, i.e. the number of electrons that pass from the cathode to the anode, is proportional to the X-ray tube current (in mAs). The application of a high-voltage yields an increase of the kinetic energy of electrons and hence their acceleration. When the electronic beam hits the anode, the electrons lose their energy through complex interactions with the nuclei of the tungsten atoms. Among all these interactions, three effects predominate: (a) heat dissipation, (b) Bremsstrahlung emission and (c) characteristic radiation.

Heat dissipation Minor interactions such as a collision with an orbital electron of the nucleus yield a transfer of energy to this secondary electron, and is eventually dissipated into heat. Studies have shown that the major part of the electron energy (at least 99%) is converted into heat. Consequently, one of the main concern in the design of the X-ray tube anode is thermal dissipation (i.e. cooling) of the anode.

Bremsstrahlung emission

The charge and mass of the nucleus are much greater than those of an electron. The corresponding energy barrier does not allow penetration in the nucleus, and the electron is consequently scattered from its initial trajectory. This deviation induces a release of kinetic energy (inelastic scattering) and the emission of a (Bremsstrahlung) X-ray photon. The energy of an X-ray photon is given by E = q • U where U denotes the voltage and q refers to the electric elementary charge: if U is set equal to 100 kV then E is equal to 100 keV, given that 1 eV is approximately equal to 1.6 • 10 -19 J. The amount of transferred energy varies such that Bremsstrahlung emission yields a continuous X-ray spectrum (in orange in Fig. 2.13), which is bounded by:

• low energy cut-off E min below which X-ray photons are so heavily attenuated that their energy can be considered negligible;

• high energy cut-off E max at the operating peak voltage. It corresponds to the maximal energy that an electron can transfer.

Characteristic radiation A collision with an electron of the inner shell of the nucleus can ionize the target atom, i.e. eject this secondary electron and leave a vacancy behind. Since this is an unstable state, the vacancy is eventually filled by an outer shell electron. This induces the emission of a quantized photon characteristic of the target material. The associated energy E C corresponds to the energy difference between the outer shell (E OS ) and the inner shell (E IS ):

E C = E OS -E IS .
These radiations have a quantization that depends on the atomic number (i.e. the chemical structure) of the target material and result in characteristic peaks in the X-ray spectrum (in blue in Fig. 2.13).

Figure 2.13: Spectrum emitted by the X-ray tube for visualization of biological tissues with an operating voltage of 150 kV and a tungsten target. The continuous spectrum of Bremsstrahlung emission is plotted in orange, while radiations characteristic of the target material are superimposed in blue.

Interactions of X-rays with Biological Tissues

At the atomic scale, interactions between X-rays and the biological tissues they penetrate mainly consist of collisions between the incident X-ray photons and the electrons surrounding the nuclei of the atoms of the tissue and result in the attenuation of the initial beam. Four independent electronic interactions exist:

• elastic scattering that occurs for very low photon energies and contributes only slightly to attenuation (less than 1%);

• photoelectric effect in which a collision with an inner shell electron of a material leaves a vacancy that may be refilled while producing a quantized X-ray photon with an energy that depends on the atomic number;

• inelastic scattering or Compton effect in which a collision with an outer shell electron of the tissue nucleus results in the transfer of a fraction of the photon energy only to the electron. A photon with the fraction of the energy retained is emitted in a deflected direction.

• pair production that is not relevant with the energies considered in X-ray CT.

Only the photoelectric and Compton effects significantly contribute to the attenuation, i.e. to the formation of the radiological contrast. The fraction of the X-ray beam that is absorbed by multiple collisions per unit thickness of material is a fundamental property of a material, that is described by linear attenuation coefficient µ. Its value basically depends on photon energy and accounts for material density and mean atomic number. In particular, biological tissues having a low atomic number emit photons of low energy, that are easily attenuated. Assuming a monoenergetic X-ray beam of energy E and intensity I 0,E , the transmitted intensity through a uniform tissue follows an exponential decrease. The measured intensity corresponds to the integration of the transmitted intensity along line D that connects the source to the measurement point:

I E = I 0,E • exp - M ∈D µ E (M )dM .
(2.12)

Typical µ E values are given in the first column of Tab. 2.1 for an energy E = 100 keV. They decrease as the X-ray energy increases, and hence the radiological contrast as well. For instance, the linear attenuation of water is found to decrease from 1.84 •10 -1 cm -1 at 80 keV to 1.51 •10 -1 cm -1 at 150 keV. The X-ray measurement is transformed into a density measurement by taking the log-transform of the ratio of the transmitted intensity given by Eq. (2.12) and the incident intensity I 0,E . The integration of the linear attenuation coefficients along D is directly related to the density measurement:

M ∈D µ E (M )dM = ln I 0,E I E . ( 2

.13)

Angiography Blood vessels have attenuation coefficients that are comparable to those of the surrounding structures. It is possible to selectively introduce contrast-enhancement agents into the vessels to modify their µ-values. This medical imaging technique is called angiography 2 and consists in injecting a radiopaque solution within the vessels, thus enhancing the blood flow. Provided that µ C refers to the enhanced linear attenuation coefficient, while µ B refers to the non-enhanced linear attenuation coefficient, the 'contrast' intensity I C that is measured after injection of the contrast agent, is simply:

I C = I 0 • exp - M ∈D µ C (M )dM = I 0 • exp - M ∈D µ C (M ) -µ B (M )dM - M ∈D µ B (M )dM = I B • exp - M ∈D vessel µ C (M ) -µ B (M )dM (2.14)
2

The word angiography comes from the Greek words angeion which means 'vessel', and graphein which means 'to write'.

where I B is the intensity measured before contrast agent is injected and D vessel determines the segment of the X-ray beam that goes through the blood vessel. Since a lower X-ray tube voltage yields stronger contrast enhancement for a given injection of contrast medium, a compromise is set between the quantity of contrast medium and the voltage to optimize the contrast enhancement.

Beam hardening

We have seen that the X-ray beam emitted by the tube is polyenergetic. The energy spectrum ω modifies Eq. (2.12):

I ω = E∈ω I 0,E • exp - M ∈D µ E (M )dM dE . (2.15)
The body tissues then act as a filter. Low energies (also called soft energies) are more attenuated than high energies (also called hard energies). This effect is known as beam hardening. Consequently, the log-transform is no longer linearly related to the linear attenuation coefficients of the encountered tissues, but instead:

ln E∈ω I 0,E dE I ω = g ω M ∈D µ ω (M )dM . (2.16)
where g ω captures the non-linearity of the measurement process and µ ω is a tissue attenuation for a given spectrum ω rather than a single energy. Note that the non-linearity is stronger if the spectrum is wider. The X-ray beam is generally filtered or 'hardened' directly out of the tube by placing a thin metal sheets (of aluminum, copper or brass), so that lower energies are preferentially filtered out and the beam becomes near monoenergetic so that remaining non-linearities can be estimated and corrected by calibrating the system.

CT for diagnostic imaging with scanners

An X-ray CT scanner such as the GE Lightspeed (GE Healthcare, Chalfont St. Giles, UK) of which a photograph is shown in Fig. 2.14 consists of a closed gantry holding the imaging system (X-ray tube, detector and high-voltage generator) that is continuously rotating around the patient, while the latter is translated perpendicularly, so that a large amount of projection data is acquired. From this acquisition, a volume of (2D) cross-sectional slices of the human body of some thickness is reconstructed. The mathematics presented in Sec. 2.1 describe the 'ideal' tomographic case, while Sec. 2.2 described the formation of the X-ray measurements. This section aims at detailing key aspects of the design of an industrial X-ray CT scanner that minimizes data corruption inherent to the X-ray physics, and collect data so that the sampling requirements are satisfied and an FBP-type reconstruction is applicable. We discuss in particular the design of key hardware components and reconstruction algorithms that are used in practice, as well as their performances.

Detector

All modern X-ray CT scanners are equipped with solid-state crystal detectors that are made from a variety of solid materials with high atomic numbers and high densities (e.g. cadmium tungstate, cesium iodide). They work on the scintillation principle, i.e. the crystal detector element fluoresces when hit by an X-ray photon and a photodiode that is attached to the crystal and affixed to a circuit board converts the photon energy into electrical energy. The key characteristics of a detector cell are:

• a high absorption efficiency (98% for solid-state crystal detectors);

• a fast decay time (30 ns for Gemstone TM the scintillator recently developed by GE Healthcare);

• a low electronic noise;

• an analog to digital conversion over 24 bits. To shorten the data collection time, all measurements of a given angle are acquired simultaneously. the detector cells with their read-out electronics are packed into a long and narrow curvilinear array. The length must accommodate all patient sizes without missing samples at a given angle (e.g. avoid truncation).

Scatter

Simultaneous acquisition results in cross-talk due to the fact that a photon that is scattered away from one detector cell is measured by another. These cross-measurements yield contrast loss. To reduce scatter and prevent contrast reduction, an anti-scatter collimator3 is placed ahead of the detector, that is 1D or 2D array of parallel lead strips that are focused towards the X-ray source. Thus, most of the scattered photons that have trajectories not belonging to a source to detector line are intercepted and absorbed by the grid. Because a part of the non-scattered X-ray photons is also intercepted by the grid, its use is accompanied with an increased dose exposure. In a complementary approach for reducing scatter radiation, as well as patient dose exposure and dynamic range, the X-ray tube of CT scanners is equipped with an additional bowtie beam shaping filter that maintains a more uniform X-ray field at the detector. It consists in compensating for the variable path length of the patient across the imaged field of view (FoV): the filter thickness increases from the center to the outer edge, which significantly reduces the dose to external parts of the patient. It results in an overall improvement in HU accuracy, low-contrast detectability and imaging dose. Different filter designs allow for imaging different regions of the body.

Dose

In practice, projection data are corrupted by quantum noise due to the limited number of X-ray photons that are recorded at the detector to form a measurement. As the flux tends toward zero, the logtransform in Eq. 2.12 becomes unreliable. Quantum noise approximately follows a Poisson distribution. Consequently, the signal-to-noise ratio (SNR) is proportional to the square root of the X-ray dose. To decrease noise, the intensity and/or the voltage of the X-ray tube must be increased, and hence the dose exposure to the patient. This is not desirable because X-rays are ionizing radiations. Although technical progress allowed for steadily reducing dose with no loss in image quality, efforts continue to further decrease the risk for the patient of developing an X-ray-induced cancer.

CT image display

The reconstructed volume is discretized with (e.g. 512×512×h where h is the number of slices considered) anisotropic volume elements (voxels) with equal, squared-shape cross-sections, that are confounded with the (e.g. 512 × 512) picture elements (pixels) of a given cross-section with submillimetric resolution. Because of the integration within the X-ray CT detector, the resolution in the longitudinal direction was historically not as good as in the transversal direction and varied from 1 to 10 mm, depending on the slice thickness, but now voxels tend to become isotropic.

Multi-planar visualization Modern scanners deliver a set of contiguous slices that form a volume that can be displayed using planes perpendicular to an arbitrary axis. There are three special axes that are illustrated in Fig. 2.15, thus resulting in axial, coronal, or sagittal visualizations of the volume.

Reconstruction is performed over planes orthogonal to the tomographic rotation axis, visualization in oblique planes is called multi-planar reformatting. Hounsfield units For a given X-ray spectrum ω, the value of the attenuation coefficient computed in tomography is quantized with value f that is expressed in Hounsfield unit (HU). By convention, the linear attenuation of a tissue is normalized relatively to the one of water and the scale is adjusted so that the value corresponding to air f air is equal to -1000 HU according to:

f (M ) = 1000 • µ ω (M ) µ ω,water -1 . (2.17)
Note that by definition f water = 0 HU. This normalization makes the reconstructed value less dependent from the energy value or spectrum, and independent from the reconstructed voxel size. An estimate of the Hounsfield units for materials of clinical interest, that was computed at 100 keV, is given in Tab. Note that with C-arm systems and in the following of this work shifted Hounsfield units (f air = 0 HU) are used by convention.

Dynamic range

The reconstructed images is displayed using a scale of 2 8 gray levels, while reconstructed values can go up to 12 or even over 16 bits in presence of metal or concentrated iodinated contrast media. Rather than scaling the original dynamic range down to 8 bits, windowing is used, that consists in interactively displaying Hounsfield units in a limited-range window that is adjusted to modify the image dynamic so that the visualization of a certain type of tissues is enhanced:

• window level (WL) sets the mean HU of the displayed window;

• window width (WW) sets the HU range coverage; Gray level G min = 0 is set to all pixels of HU below f min = W L -W W 2 while G max = 255 is set to all pixels of HU above

f max = W L + W W 2 .
A wide window is chosen to display tissues largely differing in density, while a narrow window is chosen to visualize small contrasts, as illustrated in Fig. 2.16(a) and Fig. 2.16(b) respectively.

Single-slice CT

To collect projection data, the detector of CT scanners originally consisted of a 1D array of cells arranged in a curvilinear configuration, so as to acquire a single (2D) cross-sectional slice of the human body.

Data discretization and sampling requirements

Current CT scanners use a design that appeared with the third generation of CT scanners [START_REF] Edelheit | Reconstruction of Objects from Diverging X-rays[END_REF]. It only involves the rotation of the imaging chain. At each angular position considered, a set of fan-beam lines is simultaneously collected. These are lines that diverge from a single X-ray source to cover the entire field of view of the detector. The filtered backprojection algorithm that we presented in the previous section assumes that the projection data are collected in parallel-beam geometry. Consequently, it must be adapted to handle the fan-beam geometry. To this purpose, a coordinate transformation from fan-beam geometry to parallel-beam geometry is introduced. A description of the fan-beam geometry is provided in Fig. 2.17. Let us assume the X-ray lines come from a source S at distance d SOD from the center of rotation located at the origin O of the coordinate system. Each divergent line, denoted D F B (β, γ), is determined by two parameters:

• the angle β that locates the X-ray source on its circular trajectory: β is defined as the angle between the center line of the fan and the y-axis;

• the angle γ that identifies the position of a particular ray within the fan with respect to the center line.

Let us denote p F B,β (γ) the line integral of f along divergent line D F B (β, γ), that is:

p F B,β (γ) = M ∈D F B (β,γ) f (M )dM . (2.18)
The divergent line D F B (β, γ) also belongs to the set of parallel lines, previously denoted D P B (u, φ), whose coordinates are given by the following coordinate transformation:

φ = β + γ u = d SOD • sin γ . (2.19)
The change of variable in Eq. A first reconstruction strategy consists in estimating the equivalent parallel-beam data from the data collected in fan-beam geometry in a process known as 'rebinning' [START_REF] Herman | Reconstruction from divergent beams: a comparison of algorithms with and without rebinning[END_REF]. Assuming fan-beam measurements are available at projection line D F B (β, γ), the corresponding parallel-beam measurement must be computed from the mathematical equality

p P B,φ (u) = p F B,β (γ).
Uniform discretization in (β, γ) is not uniform in (u, φ), thus requiring projection data resampling through interpolation in the (u, φ)-space, of which we do not give details here, but for which we suggest that the interested reader refers to [START_REF] Peters | Computed tomography with fan beam geometry[END_REF]. The filtered backprojection reconstruction in parallel-beam geometry is then simply applied to the rebinned data. However, because of the interpolation step, rebinning methods have to be performed on the complete projection data set, and thereby cannot be applied in a sequential manner during the acquisition, which may be computationally unattractive. In addition, linear interpolation does not fully preserve the high-frequency content of the original fan-beam projection yielding a loss in spatial resolution. In all cases, interpolation requires a large number of views, that is an angular sampling similar to the detector sampling in u because u is proportional to γ. For a limited number of views, interpolation is not adequate.

Filtered backprojection of fan-beam data

A second reconstruction strategy consists in providing an implementation of the Radon inversion formula that is appropriate for fan-beam geometry and that can be interpreted as a weighted filtered backprojection. For parallel projections, the reconstruction formula that was derived in Eq. (2.3) can be modified to include all projections over 2π:

f (x, y) = 1 2 2π 0 [p P B,φ * h](u)dφ (2.20)
Roughly speaking, corrective terms for the fan-beam geometry are introduced in the filtered backprojection algorithm [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF] by applying the change of variables from parallel-beam geometry (u, φ) into fan-beam geometry (β, γ):

f (x, y) = 1 2 2π+Γ -Γ 1 w(x, y, β) 2 • p ′ F B,β * g (γ) dβ (2.21)
This leads to the following three-step algorithm:

1. Weight the current view p F B,β (γ) with d SOD cos γ:

p ′ F B,β (γ) = (d SOD cos γ) • p F B,β (γ);
2. Filter the view by convolving the weighted view with the modified ramp kernel of impulse response:

g(γ) = γ sin γ 2 • h(γ);
3. Perform a backprojection of the filtered views, that is weighted with the squared inverse of distance w defined as the distance from the point to be reconstructed to the source

w(x, y, β) = (d SOD + x sin β -y cos β) 2 + (x cos β + y sin β) 2
We have presented in Sec. 2.1.3 a theoretical sampling requirement for tomographic systems. However, in practice, images obtained with CT scanners always present subtle subsampling streak artifacts.

Multislice CT

With a single slice CT the complete data set necessary for reconstruction of a given slice is acquired in fan-beam geometry, then the table is translated to scan the next slice and the process is repeated. Imaging of the whole body in this slice-by-slice manner requires long examination time for a patient and limited patient throughput. Increasing the volumetric coverage, i.e. acquiring more slices per second can be achieved with multislice and helical scanning approaches.

As for single-slice imaging, it is natural to collect each parallel line integral separately, this holds for three-dimensional (3D) objects that are scanned slice by slice. CT scanners allowing for the acquisition of multiple slices simultaneously (commonly 4, 8, and 16 rows) using adjacent detector arrays have been commercially available for twenty years. This scanning was developed to take full advantage of the development of new detectors made as a 2D array of detector elements, as the 3D equivalent of fan-beam scanning. The increased width of the detector allows for collecting several fan-beam projections at once, forming a cone in a so-called cone-beam (CB) geometry.

To improve the scanning time, current conventional CT scanners rely on helical-scan geometry [START_REF] Robb | Dynamic spatial reconstructor: an x-ray video fluoroscopic CT scanner for dynamic volume imaging of moving organs[END_REF]. Helical CT is a system in which the source-detectors system rotate around the patient, while the table on which the patient lies is continuously moved orthogonally to the plane of rotation. Thus the trajectory of the source relative to the patient describes a helical orbit. The principle of the helical acquisition is illustrated in Fig. 2.19. When the number of detector slices is small, all the projection lines can be assumed to be parallel to the central plane and fan-beam geometry is still relevant. Reconstruction methods rely on data rebinning for a particular reconstruction plane in the virtual equivalent projection data into parallel geometry using linear or cubic interpolation between data measured at the same angle but at different axial positions. Direct reconstruction of helical scans is beyond the scope of this presentation. The key point is that methods have been devised to adapt the FBP algorithm to helical data. Specific artifacts that are seen in helical FBP reconstruction are the so-called 'windmill' artifacts that are due to the sampling pattern.

Imaging fast moving objects

Helical CT scanners can acquire the complete projection data set related to the entire body in the range of 1-5 seconds. When imaging fast moving data, abrupt changes that may occur between the measurements result in motion artifacts as discussed in Sec. 2.1.3. An overview of the strategies that are commonly employed for reducing these artifacts is provided by Bonnet et al. [START_REF] Bonnet | Dynamic x-ray computed tomography[END_REF].

Cardiac CT

In cardiac CT, the primary goal is not to capture all beats of the beating heart, but to freeze one phase. It is thus, still, a static problem and both high temporal and spatial resolution are required, so that cardiac motion does not significantly deteriorate the image quality.

High spatial resolution of moving structures can only be achieved if data collection is fast enough to avoid motion artifact in an attempt to 'freeze the motion' during acquisition. Ultra-fast systems such as multislice CT can provide complete organ coverage in a single sub-second rotation, opening up the possibilities of functional imaging such as perfusion studies. In modern CT scanners, a 360 • acquisition of projection data can be achieved in 0.3 seconds. Short scan reconstruction only requires slightly more than 180 • . Multibeam field emission [START_REF] Lu | Multibeam field emission x-ray system with half-scan reconstruction algorithm[END_REF] have been proposed. Dual CT scanners further reduce the minimum rotation to 90 • by mounting two X-ray sources and two detectors simultaneously. The temporal resolution achieved on the dual source CT that has been recently commercialized by Siemens is claimed to be 83 ms. With technical development, the temporal resolution of multislice CT offered by manufacturers has improved from 250 ms on 4-slice scanners to 180 ms on 64-row scanners such as GE HD750. As the cone angle is increased to allow larger regions of interest to be imaged, the scatter-toprimary ratio increases significantly, which may yield image quality degradation if not handled properly. Cone-beam artifacts can be tolerated up to 16 rows, but with increasing number of rows, improved algorithms are employed, that take the CB geometry into account. To further reduce the scanning time, the width of the detector can be increased up to 128, 256 and even 320 rows. Thus, Toshiba Aquilion ONE that is equipped with a detector of 320 rows allows for dealing with imaging the entire volume in one single orbit of the CT gantry in which line integrals through all parts of the body have to be processed simultaneously.

Very fast variations within the heart lead to motion artifacts that prevent diagnostic studies and require more sophisticated clinical protocols. An essential part of cardiac imaging is to obtain a low and regular patient's heart rate to reduce the requirement for a high temporal resolution. The use of β-blockers is advocated to lower the heart rate to less than 65-70 beats per minute and to make the rhythm more regular. With advances in CT scanner technology, the heart rate range over which diagnostic studies can be obtained will likely increase. A widespread approach known as 'gating' and developed since the 1980s [START_REF][END_REF] has exploited the repeatability of the cardiac cycle to select one subset of consistent views corresponding to the same phase of the cycle and for which FBP reconstruction is feasible. To ensure that sufficient line integrals are acquired, a time window (gating window) is placed around the reconstructed phase time point. The width of this window is responsible for the temporal resolution of the system. The selection of a subset of consistent projections is done via external models that play the role of spatio-temporal indicator for the cardiac function. The electrocardiogram (ECG signal) gives information about the periodicity and offers the possibility to correlate retrospectively or prospectively -when the CT scanner acquisition is synchronized with the ECG-the views with a selected cardiac phase [START_REF] Kachelriess | ECG-correlated image reconstruction from subsecond multi-slice spiral CT scans of the heart[END_REF]. Because radiation dose is a significant issue in X-ray imaging, prospective gating is favoured rather than retrospective gating. In the prospective mode, a single phase of the cardiac cycle is imaged during multiple rotations and the data acquisition is triggered at a predetermined time window of the cycle using the ECG signal. In practice, data are frequently acquired while modulating the tube current: the peak current is used within a temporal window chosen out of the periods of significant motion, while outside this window the current is reduced, since the projection data are unlikely to be used for reconstruction.

Perfusion CT

Intravenous injection is performed to enhance the vascular structures in diagnostic imaging. After injection, the contrast medium circulates in the body, regulated by the cardiovascular system, and in particular rapidly diffuses to the brain. This results in contrast-enhanced vessels and parenchyma. Images of perfusion indexes are produced by calculating contrast uptake curves for each pixel over the whole acquisition time interval. Two problems arise [START_REF] Bae | Intravenous Contrast Medium Administration and Scan Timing at CT: Considerations and Approaches[END_REF]:

• for a short injection, the peak time-enhancement curve has a quick initial increase followed by a short peak and then quick decrease. A short injection duration (i.e. low volume or high injection rate) yields an earlier arterial peak and parenchymal enhancement. It requires a short scan delay. A fast injection would be better suited to a fast scan but requires more precise scan timing. Shortlasting contrast enhancement may be missed.

• a long scan duration requires a long contrast injection to maintain full enhancement during the acquisition: a prematurely terminated injection may yield insufficient contrast enhancement, while too long an injection may result in undesirable tissue contrast enhancement. When contrast medium is injected at a constant rate for a long injection, a typical enhancement profile for long-injection duration consists of a rapid initial rise, gradual increase, peak, and gradual decline. Patient motion may deteriorate the data. If the patient moves during the scan, an inaccurate result is produced since the contrast uptake curve no longer represents the true contrast uptake at the same location.

Motion correction is relatively simple in the case of brain perfusion, since a rigid registration can be performed over all reconstructed images to ensure proper location registration. The problem is much more complex when dealing with body organ perfusion since the motion is non-rigid. Perfusion over large volumes such as the liver is therefore a work-in-progress. Fast scan speed and coverage with multislice CT allow for a complete acquisition in the clinical routine in less than 10 seconds and facilitate acquisition at multiple precisely defined phases of contrast enhancement and may improve the detection and characterization of lesions in various organs. Optimization of phase-specific contrast enhancement in an organ is however complex and would require higher performance of the CT scanners than the current state-of-the-art.

Software temporal management

The technological progress to speed up the data acquisition of CT scanners or increase its coverage allow for neglecting the motion effects in most of the diagnostic applications of CT scanners. However, in spite of these advances in fast CT, it is impossible to remove all motion artifacts: among others, it has been shown that scanning time less than 50 ms are not adequate to resolve respiratory-induced motions of the pulmonary vessels [START_REF] Ritchie | Minimum scan speeds for suppression of motion artifacts in CT[END_REF] and results in intra-slice inconsistencies. So far, we considered CT hardware improvements only. Software solutions have received considerable attention in the academic research. Motion-compensated reconstruction methods mostly rely on an initial subsampled motionfree reconstruction. The analysis and extraction of parameters of cardiac shape and function allows for exploiting all the projective measurements so that one may achieve both better SNR and reduced motion artifacts. In this way, the advantages of both gated and ungated methods are kept. Common strategies rely on modeling the motion (e.g. with parametric models or dense vector fields) so that each phase can be deduced from a reference one by the motion model. A general review of 3D modeling for functional analysis of cardiac images in different modalities is given by Frangi et al. [START_REF] Frangi | Three-dimensional modeling for functional analysis of cardiac images, a review[END_REF]. Motion estimation is a challenging estimation problem though. Due to their high complexity, the introduction of the motion models in the reconstruction is not a trivial task and is still an active area of research. It has been shown by Desbat et al. [START_REF] Desbat | Compensation of some time dependent deformations in tomography[END_REF] that deformations that preserve the acquisition line geometry, i.e. maintain the acquisition X-rays as straight lines, can be analytically compensated by redefining a virtual trajectory of the source. However, elastic deformations associated to the cardiac motion transform straight lines into curves and thereby do not satisfy the linearity preservation requirement. In another work from Taguchi et al. [START_REF] Taguchi | Motion compensated fan-beam reconstruction for nonrigid transformation[END_REF], an approximate algorithm was shown to provide significant reduction of the motion artifacts due to non-rigid motions within the imaged object, but could not deal with the arbitrary non-rigid deformations associated to the cardiac motion. However, due to the complexity of the cardiac motion, 3D model-based approaches lack of robustness and rigorous clinical validation and until now very few industrial applications have followed 4 .

Chapter 3

X-ray interventional imaging

Medical imaging has known great technological advances over the past decades to become a powerful tool for the current clinical practice. On one hand, diagnostic imaging techniques have facilitated the detection, characterization and follow-up of many pathologies. On the other hand, interventional radiology (IR) that involves minimally-invasive clinical procedures has developed tremendously. In these procedures, medical devices are inserted and manipulated under image guidance through the vascular system to the pathology location, where they are used to deliver the treatment. The very first imageguided intervention (a percutaneous vascular catheterization) was proposed by Seldinger in 1952 [START_REF] Seldinger | Catheter replacement of the needle in the percutaneous arteriography: A new technique[END_REF].

Contributions to the field of IR include the works of Dotter (percutaneous revascularization, 1964) [START_REF] Dotter | Transluminal treatment of arteriosclerotic obstruction[END_REF], Grüntzig (percutaneous angioplasty, 1977) [START_REF] Grüntzig | Non-operative dilation of coronary artery stenosis. Percutaneous transluminal coronary angioplasty[END_REF] and Cato (chemoembolisation, 1981). The use of IR as therapeutic care for cardiac, vascular or oncological purposes has never stopped from increasing since then. C-arm systems, that are designed for real-time guidance (2D projective images with low-dose X-rays, also known as fluoroscopy) of interventional procedures, deliver tomographic images as well. The introduction of 3D visualization of the vascular anatomy [5] and adjacent structures through X-ray CT for interventional procedures has allowed for evaluating ambiguous anatomies and position accuracy of interventional tools, and even detecting haemorrhage. The specificities of C-arm systems with respect to X-ray CT scanners are discussed in Sec. 3.1. In order to better understand the sampling requirements that are associated to interventional data, this chapter gives a rapid overview of current angiographic procedures in Sec. 3.2; the focus is then put on embolisation procedures that consist in inducting vascular occlusion by introducing an embolic agent into a vessel through a selectively placed catheter for therapeutic purposes in the brain and in the liver in Sec. 3.3. For both cases, the temporal variations as well as the resulting image quality are discussed in Sec. 3.4.

CT for interventional imaging with C-arm systems

Although closed gantry CT scanners are the modality of choice for diagnostic volumetric imaging [START_REF] Kalender | X-ray computed tomography (review)[END_REF], they are only rarely used during minimally-invasive interventional treatments due to the ring structure which restricts the access to the patient during the procedure. C-arm systems on the other hand were developed to provide real-time 2D X-ray images of the patient for procedure guidance and hence designed as an open structure to give the medical team full access to the patient. A C-arm unit is designed such that all components within the system are carefully optimized to provide an image quality matching the clinical needs of an intervention and improve the ease of use. The key components are:

• the mechanical gantry and the table on which the patient lies;

• the X-ray imaging chain;

• the monitors that are carried by telescopic arms and used to display the images in real-time;

• and the workstation that is used to store and manipulate the images.

There is no one-fit-for-all system and C-arms come in many flavours. Mobile C-arms (see Fig. 3.1(a)) are the smallest systems and are mostly used in invasive surgery at specific steps of a procedure and moved away otherwise. They can be shared between several operating rooms. The counterpart to mobility is limited capabilities; in particular, they rarely feature computed tomography. For minimally invasive vascular procedures, fixed-room systems are used, such as the GE Innova 4100 which was used to acquire some of the data presented in this manuscript and of which a photograph is given in Fig. 3.2. It is a floormounted system. Biplane systems (see Fig. 3.1(b)), such as the IGS 630, provide two simultaneous image chains by adding a ceiling-mounted system to the floor-mounted one. Ceiling C-arms are mounted on rails to translate them in and out of the operating field. The latest designs combine features and image quality of a fixed room with increased mobility. For instance, the Siemens Artis Zeego (see Fig. 3.1(c)) is a robotmounted C-arm which is fixed and floor mounted, but powerful enough to move the C-arm away when needed. The GE IGS 730 (see Fig. 3.1(d)) is also a robot-mounted C-arm, where the robot is equipped with wheels and laser guidance so as to provide the mobility of a mobile unit within the operating room with millimetre precision in the repositioning. On 'fixed' systems, computed tomography has become essential, but still comes second in terms of requirements with respect to fluoroscopy. This section thus aims at highlighting the key differences that make C-arm systems different and less performing than CT scanners for tomographic reconstruction.

Characteristics of the C-arm acquisition

For interventional procedure guidance, interventional imaging involves several intraarterial injections that enhance the blood vessels and allow the clinician to manipulate the interventional tools through the patient's arterial system.

Mechanical considerations

The design of an angiographic system aims at fullfilling the specificities of the interventional context: (a) an easy positioning with space to work around and (b) a large freedom of movement, allowing acquisition of X-ray images (views) with various angles. The products manufactured by GE Healthcare thus benefits from various degrees of freedom. The imaging system can be oriented along three degrees of freedom through a LC architecture illustrated in Fig. 3.2:

• The L-arm holds the C-arm and can be rotated along one (L) axis, while the C-arm that holds the X-ray tube and detector can be rotated about two axes:

-the CRAnial/CAUdal angle that describes rotation of the C-arm in the patient's sagittal plane;

-the LAO/RAO: the Left (respectively Right) Anterior Oblique angle that describes rotation of the C-arm in the patient's axial plane;

• Fixed-room C-arms are isocentered: all rotations are made with respect to a common point, the isocenter O. Therefore the table is made with three additional degrees of freedom in translation so that the anatomy of interest can be put into the X-ray beam;

• Finally, the detector is equipped with a lift that modifies the source to detector distance d SID , which changes the image-chain magnification, or equivalently, the size of the field of view.

Flat-panel detector

Three types of X-ray detectors have been used in angiography systems:

• phosphor and fluorescent screens that were historically used by the very first angiography systems 1 ;

1 The term 'fluoroscopy' that denotes the resulting 2D X-ray images comes from there.

• image intensifiers;

• and digital flat panels that equip current products such as the IGS 540.

Until ten years ago, angiography systems relied on an image intensifier which consists of an electronic device that converts X-ray photons into electrons to form the digital image by photoelectric effect. Its introduction provided a far better dynamic range than fluorescent screens. Still, images suffered from distortions, and the field of view was limited to a disk. The progress made in the development of sensor systems lead to the replacement of image intensifiers by digital flat-panel detectors that provide distortion-free images, higher dynamic range and enable fast digital read-out with a maximum of 50 fps for IGS systems. Thanks to their higher DQE (detector quantum efficiency), digital flat-panels enable dose reduction or image quality improvements. The size of GE detectors depends on the imaged anatomy:

• in cardiac imaging a small 20 cm × 20 cm square panel is enough;

• biplane systems for neuroradiology use two medium 30 cm × 30 cm square panels;

• while for peripheral imaging (extremities, abdomen) a large 40 cm×40 cm square panel is preferred.

Rectangular panels of 30 cm × 40 cm are proposed instead by Siemens or Philips. The working principle of GE's detector is schematically illustrated in Fig. 3.3(a). X-ray photons are converted into visible photons by a crystalline cesium iodide (CsI) scintillator. These visible photons are in turn converted into electrons using an amorphous silicon (a-Si) photodiode to form a digital image by photoelectric effect.

Collimation

Patient radiation exposure concerns are crucial in X-ray imaging. In interventional procedures, it is also crucial to monitor clinician scatter exposure. Protective rules enable lower scatter radiation exposure: the clinician wears protective clothing and ideally supervises the scanning procedure from the control room or at least stands behind a protection glass panel. Restraining the FoV that is provided by the source-detector system so that only a region of interest (RoI) is imaged enables both dose and scatter reduction. Thereby, it seems very reasonable to limit the dose exposure to the anatomical volume that need to be imaged only. A collimator is used to restrict the X-ray beam to a specific FoV: blades that are opaque to radiation are placed on top of the X-ray tube. They can be translated to frame the target region, the complementary regions being shielded by the blades. To further reduce scatter in the measurements, a removable 2D grid is placed ahead of the flat detector.

Tomography with C-arm systems

A sequence of views (scan or spin) can be generated by a single rotation of the C-arm in the axial plane, where the angle LAO/RAO varies from -100 • to +100 • , while CRA/CAU and SID are kept constant (CRA/CAU is generally set equal to zero). The rotation over 200 • is called short-scan by opposition to a full 360 • rotation as in conventional CT. Figure 3.4 shows a few views from a rotational scan: each view corresponds to a particular position of the C-arm imaging system and, thus, illustrates a particular (2D) point of view on the (3D) imaged object.

Rotation speed and frame rate

The frame rate of the system depends on the design of the read-out electronic of the detector: it is limited to 30 fps with Innova systems and to 50 fps with IGS systems. The choice of the rotation speed depends on the clinical protocol. The C-arm gantry is used at maximal rotation speed (40 /s with IGS systems). A reminder of these rotation speeds is given in Tab.3.1.

Cone-beam X-ray transform

With a square flat-panel detector, the projection lines belong to a cone. The CB projection transform p CB,β (u ′ , v ′ ) associates to any point M of the imaged FoV the integral projection of function f along line D CB (u ′ , v ′ , β) as illustrated in Fig. 3.5:

p CB,β (u ′ , v ′ ) = M ∈DCB(u ′ ,v ′ ,β) f (M )dM . ( 3

Projection matrix

Projection matrices model the set of line integrals defined by the geometry of acquisition. It is a fourdimensional (4D) space, while the set associated to the detector plane is a 3D space. Thus, the projection transform 3DX is described via a 3 × 4 projection matrix [START_REF] Faugeras | Three Dimensional Computer Vision: A Geometric View-point[END_REF] considering three different coordinate systems: (a) the real world referential, (b) the X-ray source referential, and (c) the detector referential. M , y ′ M , z ′ M ) the absolute position of point M of the imaged FoV, in the 3D world coordinate system. Let us assume M of coordinates (x M , y M , z M ) in the X-ray source referential is projected on the detector at the position in the (2D) detector referential that is given by the coordinates (u ′ , v ′ ).

Intrinsic parameters

The intrinsic parameters characterize the projective X-ray process. They are independent from the X-ray source point of view and encode for:

• the magnification induced by the projection process, that is associated with the relative position of the source with respect to the detector α = d SID ;

• the detector orientation which corresponds to the (orthogonal) projection of the X-ray source (u S , v S ).

Detector

Optical axis

α = d SID X-ray source S (u S , v S ) u ′ v ′ Real world referential O x y z x ′ y ′ z ′ M Figure 3
.6: Cone-beam geometry. Ideal acquisition: x, z, x ′ , y ′ and u ′ belong to the same plane and z ′ , v ′ and y are aligned.

The coordinates of the projection of M on the detector can be written as a linear mapping between homogeneous coordinates:

  su ′ sv ′ s   =   α 0 u S 0 α v S 0 0 1   •   x M y M z M   (3.2)
From this equation, we observe that the magnification factor is modulated by coordinate s = z M that characterizes point M location along the optical axis so that local magnification takes value α s .

Extrinsic parameters

The extrinsic parameters identify the transformation between the X-ray source referential and the real world referential. Using the extrinsic source parameters, it is straightforward to find the relation between the absolute coordinates of point M and its relative coordinates in the X-ray source referential. The rigid transformation between the two referentials is the combination of:

• a 3D rotation matrix that encodes for the orientation of the X-ray source;

• a 3D translation vector that encodes for the position of the source within the 3D real world coordinate system. The origin of the system is set arbitrarily. For CT, it is fixed at the isocenter (i.e. the center of rotation).

In the ideal case, the trajectory of the C-arm is circular (i.e. the source rotates around the z-axis only). Hence, the coordinates of the points that belong to line D CB (u ′ , v ′ , β) are defined by:

  x M y M z M   =   cos β sin β 0 0 0 0 -1 0 -sin β cos β 0 d SOD   •     x ′ M y ′ M z ′ M 1     (3.3)
With these notations, the projection matrix is the product of the matrix of the intrinsic parameters that is given in Eq. 3.2 with the matrix of the extrinsic parameters that is given in Eq. 3.3. The projection operator is completely described by the set of projection matrices, where one matrix encodes for the geometry of a view (i.e. for a given angle):

  su ′ sv ′ s   =  
  •     x ′ M y ′ M z ′ M 1     . (3.4)
Since the trajectory of the C-arm deviates from the ideal circular orbit, coefficients p ij are directly estimated. The computation of a CB reprojection or backprojection with these matrices is described in [START_REF] Riddell | Rectification for cone-beam projection and backprojection[END_REF].

Cone-beam reconstruction

One of the main issue of CB reconstruction with C-arm is that the source trajectory is limited to a circle. This configuration does not allow for exact reconstruction, meaning that there will always be an error independently from the amount of views and from the detector resolution. Tuy [START_REF] Tuy | An Inversion Formula for Cone-Beam Reconstruction[END_REF] proved that exact reconstruction is feasible in cone-beam geometry for slices that are crossed by the trajectory of the X-ray source. Smith [START_REF] Smith | Image reconstruction from cone-beam projections: Necessary and sufficient conditions and reconstruction methods[END_REF] then proved that this condition is not only sufficient, but necessary. In 1984, Feldkamp, Davis and Kress published the so-called Feldkamp algorithm (FDK) method [START_REF] Feldkamp | Practical Cone-Beam Algorithm[END_REF], a practical FBP algorithm for cone-beam geometry. Although many variants and alternatives [START_REF] Herman | Mathematical framework of cone beam 3d reconstruction via the first derivative of the radon transform[END_REF] to the original FDK method have been suggested, FDK algorithm is still the standard reconstruction for all C-arm systems. In comparison to a regular FBP implementation, FDK requires only slight modifications of the filtering and backprojection steps. In the central slice, the FDK algorithm reduces to FBP for the fan-beam geometry and therefore FDK produces exact reconstruction in this transaxial slice, that is the only one crossed by the source trajectory.

Fan-beam reconstruction for flat detectors

The reconstruction results of fan-beam geometry with a flat detector is a preliminary task for understanding the Feldkamp algorithm. In Sec. 2.3.2 we assumed the projection data to be sampled in equiangular intervals on a circular detector. However, C-arm systems use equidistant pixels on a flat detector. Similar to the case of equiangular fan beams, a coordinate transformation can be applied to the FBP integral between the fan-beam parameters (β, u ′ ), where u ′ refers to the abscissa of the projection of M on the flat detector, and their parallel-beam equivalent (φ, u). As seen in Fig. 3.7, the relation is:

         φ = β + γ = β + tan -1 u ′ d SID u = u ′ cos γ = u ′ d SID d 2 SID + u ′ 2 ,
The resulting FBP algorithm then becomes:

f (x, y) = 1 2 2π 0 1 w(x, y, β) 2 • p ′ FB,β * h (u ′ ). dβ. (3.5)
This leads to the following three-step algorithm:

1. Weight projection p FB,β (u ′ ) according to the magnification factor cos γ = d SID

d 2 SID + u ′ 2 : p ′ FB,β (u ′ ) = cos γ • p FB,β (u ′ ); 2.
Filter the projection by convolving the weighted projection with the ramp kernel;

3. Perform a backprojection of the filtered projection, that is weighted with the ratio between projection onto the optical axis and the source center distance d SID :

w(x, y, β) = (d SID + x sin β -y cos β) d SID .

Cone-beam reconstruction for flat detectors

Cone-beam reconstruction is carried out through the FDK algorithm that is defined as:

f (x, y, z) = 1 2 2π 0 1 w(x, y, β) 2 • p ′ CB,β * g (u ′ , v ′ ). dβ (3.6)
Like the reconstruction approach in fan-beam geometry, the FDK algorithm requires three successive steps:

1. Weight a projection according to its position within the 3D cone, that corresponds to the cosine of the angle between the considered projection measurement and the central ray of the current projection:

p ′ CB,β (u ′ , v ′ ) =   d SID d 2 SID + u ′ 2 + v ′ 2   • p CB,β (u ′ , v ′ );
2. Filter the weighted view along the horizontal detector lines with the ramp kernel;

3. Perform a weighted backprojection of the filtered projections along the cone with the same weighting factor w(x, y, β) that was used in the 2D fan beam case.

Note that similarly to the reconstruction in fan-beam geometry (see Sec. 2.3.2.1), instead of using FDK algorithm, a rebinning approach adapted to cone-beam reconstruction can be used. A set of virtual parallel X-ray beams with different angles for different detector rows is generated from rebinning of cone-beam data.

Short-scan algorithm

With C-arm systems, it is not possible to rotate the X-ray tube and the detector around the patient by a full 360 • . A partial rotation over 180 • plus twice the fan angle Γ provides sufficient information for reconstruction of the data but also introduces some redundancy. As shown in Sec. 2.3.2.1, the reconstruction algorithm has to be adapted to incorporate the short-scan geometry. A classic design consists in zeroing the data in a triangle of the (β, γ)-diagram plotted in Fig. 2.18, but it introduces sharp edges between non-redundant and redundant data. These edges are enhanced by the filtering with the ramp filter and lead to strong artifacts in the reconstruction. To smooth this discontinuity, Parker's weights [START_REF] Parker | Optimal short scan convolution reconstruction for fan-beam CT[END_REF] that rely on a Parker window function as defined in Eq. 3.7 are used:

w β (γ) =            sin 2 π 2 • β Γ -γ if 0 ≤ β ≤ Γ -2γ 1 if Γ -2γ ≤ β ≤ π -2γ sin 2 π 2 • 2π + 2Γ -β Γ + γ if π -2γ ≤ β < π + Γ . (3.7)
The Parker's weights are applied to an entire column prior to ramp filtering to incorporate the short-scan geometry into the FDK reconstruction.

C-arm CT image quality

Image quality can be defined in terms of artifacts, noise, temporal and spatial resolution, as well as contrast resolution. However, image quality is not determined by the features of the detector array only, but also by the nature of the imaged data and the variables selected by the operator. The image quality is qualitatively assessed by visual inspection focusing on both the depiction of anatomic structures of clinical interest and the introduction of artifacts by the reconstruction process. A summary of the performance of current C-arm systems for the tomographic reconstruction is based on the examination of its most significant parameters is given in Tab. 3.1.

Artifacts

We have introduced in Sec. 2.1.3 the subsampling and motion streak artifacts. Because of the wide variety of artifacts that may be encountered in X-ray CT imaging, there exists no standard definition of what is an artifact. However, for the needs of this manuscript, we will define an artifact as any content that is created by the reconstruction process and that does not represent the distribution of the object being imaged except noise, i.e. any defect in the image that cannot be removed by increasing the dose. Artifacts originate from insufficient correction of the physical effects that we just described and sum up with the sampling artifacts described in Sec. 2.1.3. It is thus not always obvious to identify the root cause of an artifact and one issue might be actually hidden by another, so that correcting for one artifact only reveals another one. In the following we shall only address correction of subsampling artifacts. Because the maximum rotation speed (40 • per second) will never meet that of a CT scanner (superior to 1000 • per second), the temporal resolution is not better than 4 seconds. Image quality is thus degraded by motion artifacts as we shall discuss in Sec. 3.4.2.

Low contrast visualization

Contrast resolution is affected by the higher scatter-to-primary ratio and by the sampling limitations.

As the cone angle is increased, the scatter-to-primary ratio increases: thus, the ratio goes from 14% for a cone angle of about 0.5 • to greater than 120% for a cone angle of 7 • [127]. This is illustrated in Fig. 3.8. The limitations of the anti-scatter device of C-arm systems (see Sec. 3.1.1.2) leave a significant contribution of the scatter within the data [START_REF] Gupta | Ultra-high resolution flat-panel volume ct: fundamental principles, design architecture, and system characterization[END_REF]. Efficient scatter suppression and additional correction procedures are essential for C-arm CT to achieve CT-like image quality. Scatter correction approaches include measurement techniques that require additional hardware, software models that operate directly on the projection data [START_REF] Zellerhoff | Low contrast 3d-reconstruction from c-arm data[END_REF], and hybrid approaches [START_REF] Siewerdsen | A simple, direct method for X-ray scatter estimation and correction in digital radiography and cone-beam CT[END_REF][START_REF] Zhu | An efficient estimation method for reducing the axial intensity drop in circular cone-beam ct[END_REF]. Iterative approaches appear promising to obtain further improvements beyond state-of-the art scatter correction methods [START_REF] Kyriakou | A simple, direct method for X-ray scatter estimation and correction in digital radiography and cone-beam CT[END_REF]. In addition, subsampling due to limited framerate and CB circular trajectory [3] degrade low contrast detection through the mecanism illustrated in Sec. 2.1.3. Another issue is data truncation. The maximum size of a flat panel detector is 40 cm. Given a system magnification of about 1.65, the maximum FoV is limited to 24 cm. Thus, in most of the examinations, the object of interest is only partially exposed and is said to be truncated. The reconstructed RoI is composed of voxels that are observed by each source position.

Truncation introduces a strong cupping that is compensated by data extrapolation at truncated edges. This correction leaves a residual cupping that alters the visualization of small contrast as illustrated in Fig. 3.9. To summarize, strategies to optimize low contrast detection involve:

• small vertical cone angle to decrease scatter;

• scatter subtraction after measurement under the blades;

• slow rotation with large horizontal cone angle that enables avoiding truncation to decrease subsampling artifacts.

A comparative evaluation of C-arm CT image quality that is shown in Fig. 3.10(a) with respect to image quality with CT scanners that is shown in Fig. 3.10(b) emphasizes that differentiation of the brain soft tissues is easier due to higher contrast resolution. The contrast resolution for Innova is evaluated in Fig. 3.10(c).

Vessel visualization

Vessel visualization requires a high spatial resolution, that can be much better with C-arm systems (voxel size of about 0.2 mm) than with CT scanner (voxel size up to 0.6 mm). However, the flat panel scintillator suffers from greater lag (afterglow), which results in a slower response than the CT ceramic detectors [START_REF] Gupta | Ultra-high resolution flat-panel volume ct: fundamental principles, design architecture, and system characterization[END_REF] and limits the frame rate. Special read-out techniques such as pixel binning can be applied to increase the speed. It consists in combining the outputs of neighboring pixels of the detector into one single reading: for instance, 2 × 2 binning implies that a total of 2 × 2 neighboring pixels are read together. In return, spatial resolution is reduced since pixel binning modifies the effective pixel size.

Taking full advantage of the high spatial resolution of digital flat panels would lead to low detector framerate, which would restrict its clinical use to some high-contrast applications only. C-arm systems equipped with flat-panel detectors provide better spatial resolution but still lower dose efficiency (and thus lower contrast resolution) when compared to conventional CT systems.

To summarize, strategies to image vessels involve:

• intraarterial injection providing high contrast;

• fast rotation to limit injection time;

• small (truncated) FoV to avoid pixel binning.

Cross-sections are not very informative for the visualization of tubular structures such as the blood vessels. Maximum intensity projection (MIP) images that were first developed for use in nuclear medicine by Wallis et al. [START_REF] Wallis | Three-dimensional display in nuclear medicine[END_REF], are produced by simulating the projection -generally in parallel-beam geometry-of the voxel that has the maximum intensity only. By changing the angle of the simulated beam, a series of 2D MIP images can be generated. Maximum intensity projection is particularly adapted to the visualization of the enhanced vessels and has been shown to be quite robust to subsampling artifacts that may deteriorate the original volume. An example of MIP representation is provided in the second row of Fig. 3.13. Another suitable alternative to cross-section visualization is volume rendering (VR). It is a 3D rendering of the contrast-enhanced vessels only. It is the visualization of choice for faithfully representing small vessel abnormalities such as aneurysms [5] as can be seen in Fig. 3.14.

Obviously, it is currently impossible to provide low contrast detection of the soft tissues while injecting vessels with a single acquisition on a C-arm system.

Reference for visual comparison of clinical data

In numerical simulations, we can judge the image quality by comparing it with a ground-truth that is provided by the numerical phantom. When dealing with clinical data, there exists no ground-truth. In this work, a reference reconstruction technique is chosen, from which establish the advantages of our algorithms. Note, however, that there are other approaches to determine the quality of a reconstruction algorithm, e.g. task-based metrics [START_REF] Gang | Analysis of Fourierdomain task-based detectability index in tomosynthesis and cone-beam CT in relation to human observer performance[END_REF] or visual quality assessment by experts.

Normalized root mean square deviation The convergence behaviour and quality of the final solution were monitored in numerical simulations through the calculation of a similarity measure for each reconstruction result. In this work, we report on the normalized root mean square deviation (RMSD) measure d r between a given reconstruction result f and the ground-truth f over several RoIs of K r pixels in the background structures, that emphasizes the importance of a few large errors. It is mathematically defined as:

d r = Kr k=1 (f ) k -( f ) k 2 Kr k=1 ( f ) k 2 .
RMSD may not be the best indication for image quality. Note that other measures exist, although not used in this work. They emphasize other aspects of image quality, such as:

• the normalized mean absolute deviation, that emphasizes the importance of many small errors;

• or the structural similarity index (SSIM), that emphasizes the degradation of structural information [START_REF] Wang | Image Quality Assessment: From Error Visibility to Structural Similarity[END_REF].

Angiographic procedures

The standard clinical interventional workflow that is summarized in Fig. 3.11 involves (a) a pre-procedural evaluation step that relies on diagnostic imaging in order to decide and plan the intervention, (b) the intervention by itself that is performed under local or general anaesthesia and X-ray C-arm guidance, and (c) a post-procedural diagnostic imaging step in order to evaluate the success of the intervention.

Pre-procedural evaluation

Arterial angiography

Intravenous injection of a contrast medium can be carried out by trained dedicated paramedical personnel. Intraarterial injection, such as in the femoral artery, presents a risk of haemorrhage and can only be performed by interventional radiologists, which is today a dedicated medical speciality. A system of guidewires and catheters are introduced within the arterial system and vessel enhancement is obtained through the injection of iodinated contrast agents. This technique allows for local enhancement of vessels with higher contrast than intravenous injection.

To operate within the arterial system, access to the blood vessels is performed most frequently through the femoral artery. The access technique begins with placing a radiopaque marker (such as a clamp) at the point of maximal arterial pulse and making a small incision through which a needle is then introduced. Once the needle tip is entered into the artery, a guidewire is slightly advanced. Fluoroscopic guidance ensures the proper positioning of the guidewire. The needle is then removed and the interventional tools (such as a catheter which is a long, thin, flexible and hollow plastic tube) can be introduced safely over the guidewire. The guidewire is pushed through the arterial system under fluoroscopic guidance until the pathology location is reached. Then, the interventional tools are delivered over the wire.

Digital subtraction angiography

Depending on the injection point and the amount of dye, selective parts of the vasculature can be imaged to provide to the physician a 'road-map' showing the guidewire to help guidance. However, it might be difficult to separate the vascular structures from surrounding bones or dense devices such as coils. To suppress the background structures, views are acquired under the same angle using the principle of digital subtraction angiography (DSA) [START_REF] Brody | Digital subtraction angiography[END_REF]. The first view, called mask and shown in Fig. 3.12(a), is acquired without injection, while all other views, called contrast and of which an example is given in Fig. 3.12(b), are acquired after vessel contrast enhancement. All structures but the vessels are removed by digital log-subtraction of the mask from the contrast views as shown in Fig. 3.12(c). A major cause of image quality degradation in subtracted protocols is the motion that may occur between acquisition of the mask and the contrast scans. With the progresses of X-ray flat panels, tomographic acquisition (scan) of 2D X-ray projection views has become standard to reconstruct a 3D model of the injected vessels in their environment. The principles of DSA also applies in 3D: a mask scan and a contrast scan are acquired in a single protocol. From these scans, one can reconstruct three different volumes consisting of:

• the mask volume of bones, tissues and devices that represents the vessel context alone as illustrated in Fig. 3.13(a);

• the contrast volume that displays the vessels and their context as illustrated in Fig. 3.13(b);

• the subtracted volume that displays the vessels alone as illustrated in Fig. 3.13(c).

To reconstruct a satisfying subtracted volume, the mask and contrast scans are acquired with identical parameters, which allows the straightforward removal of the redundant background structures within the projections so that they do not interfere with the reconstruction of the vessels.

Clinical use of rotational angiography for embolization

The focus of this thesis is put on embolisation procedures for therapeutic purposes in the brain and in the liver.

Cerebral Aneurysm embolisation

Rotational angiography is particularly useful for visualizing vessel stenoses, aneurysms and arteriovenous malformations (AVM)2 . Here we look at the specificities of aneurysm imaging in order to highlight the impact of sampling. An aneurysm is a vascular pathology in which an abnormal dilation develops in an artery because the wall of the vessel is weak: a normal artery wall is made up of three layers, while the structural support of the muscular layer is absent in a pathological artery. This most often results in a saccular aneurysm as illustrated in Fig. 3.14. When located in the brain, a widening aneurysm may compress the surrounding nerves and brain tissues, accompanied with different symptoms such as nerve paralysis, severe headache, pain as well as nausea [START_REF] Brisman | Cerebral Aneurysms[END_REF]. Note however that 90% of the cerebral aneurysms are present without any symptoms. The rupture of a cerebral aneurysm is usually sudden. The resulting bleeding may damage To decide on whether or not the aneurysm must be treated, the risk of rupture must be evaluated.

The growth and rupture of an aneurysm can be explained mechanically and pathologically on the basis of wall structure and blood flow. However, the most significant predictors of the risk of rupture are still unclear. For unruptured cerebral aneurysm, the predictors that are used in the clinical practice are the size and location of the aneurysm. Several studies [START_REF] Wiebers | Unruptured intracranial aneurysms: natural history, clinical outcome, and risks of surgical and endovascular treatment[END_REF] have shown that smaller aneurysms present a lower risk of rupture. Particular locations in the brain have been shown to be more likely to rupture, regardless of the size of the aneurysm. Note that around 90% of the cerebral aneurysms are less than 10 mm in diameter. Endovascular coils can be used to occlude the aneurysm so that the blood flow is no longer directed to the aneurysm. This technique is an alternative to surgical clipping, that is associated with a shorter recovery time. It is adapted to aneurysms that present a 'narrow neck'. A catheter is advanced to the location of the aneurysm. Once the catheter is inside the aneurysm, detachable small, soft metal (generally platinium) coils [START_REF] Guglielmi | Electrothrombosis of saccular aneurysms via endovascular approach. Part 1: Electrochemical basis, technique, and experimental results[END_REF] are inserted through the catheter to fill the aneurysm sac. Note that the coil size has to be adapted to the size of the aneurysm sac. The coils initiate a healing reaction, also known as clotting or thrombotic reaction, that helps block the flow of blood into the aneurysm, preventing it from further rupturing or leaking. Depending on the morphology, a stent can be deployed in order to prevent coil herniation3 . Aneurysm imaging requires a high spatial resolution as shown in Fig. 3.14 where it is fundamental to visualize the two vessels that must not get occluded by the treatment. Contrast dilution within large aneurysm is incomplete and yields artifacts. Post-procedure imaging must assess that blood does not penetrate the metal-filled aneurysm any more. In cerebral imaging, patient motion is either negligible if the acquisition is made short enough or may be compensated through a rigid registration. Consequently, cerebral vessel imaging was the first application of DSRA. In particular, several works [START_REF] Bidaut | Second-generation three-dimensional reconstruction for rotational three-dimensional angiography[END_REF][START_REF] Heautot | Analysis of cerebrovascular diseases by a new 3-dimensional computerised X-ray angiography system[END_REF]4] demonstrated the feasibility of DSRA for endovascular treatment of intracranial aneurysms, in particular for planning the treatment, i.e. finding a working view, and performing accurate measurements, as well as its superiority with respect to DSA. In particular, 'roadmapping' that consists in combining the 3D vasculature of Fig. 3.16(b) with the fluoroscopy of Fig. 3.16(a) and resulting in the merged image presented in Fig. 3.16(c) facilitates guidance of the procedure.

Transarterial chemoembolisation (TACE)

Hepatocellular carcinoma (HCC) is the most common liver cancer and has a poor prognosis. Surgical resection is the treatment of choice for patients with well-preserved hepatic function. For the treatment of unresectable tumours, transarterial chemoembolisation can be used. It is a procedure that involves intraarterial chemotherapeutic drug delivery directly into a tumour through its feeding blood vessels.

The intervention causes the tumour necrosis but also has a synergistic effect because embolic agents that are injected via the catheter subsequently block the supplying vessel (ischemia), thus decreasing the tumour supply in nutriments. The key advantages of transarterial chemoembolisation compared to systemic chemotherapies is the high local concentration (pharmacological studies indicate a 10-to 50-fold increase) of chemotherapeutic drug to be in contact with the tumour for a prolonged time, that maximizes the cell necrosis while minimizing the systemic toxicity. The efficacy of transarterial chemoembolisation to treat liver tumours was assessed by two randomized trials [START_REF] Rougier | Efficacy of chemoembolization for hepatocellular carcinomas: experience from the Gustave Roussy Institute and the Bicêtre Hospital[END_REF]. The normal liver parenchyma gets a dual blood supply from

• the portal vein that contributes up to 75-80% of hepatic blood flow,

• the hepatic artery that contributes the remaining 20-25% [START_REF] Sahani | Imaging the liver[END_REF][START_REF] Namasivayam | Hypervascular Hepatic Focal Lesions: Spectrum of Imaging Features[END_REF],

Some tumours such as HCC receive their blood supply primarily from the hepatic artery, with limited supply from the portal vein. Consequently, hepatic artery embolisation causes more ischemia to the liver tumour while the remaining normal hepatic parenchyma obtains sufficient oxygenation from the portal venous blood supply to preserve the liver function.

The standard hepatic vascular anatomy is quite complex as can be seen in the schematic representation given in Fig. 3.18. It is primordial to identify the main vessels that specifically supply the tumour, that are referred to as 'feeder vessels'. This is however not an easy task because:

• there is a large variability of the hepatic vascular anatomies from one individual to another: a patient may have only one, or two and even three hepatic arteries and the vascular territory of hepatic arteries can vary from small (subsegmental) to large (lobar) distribution, which makes difficult the accurate segmental localization of the tumour;

• it is essential to determine whether a collateral blood supply is present, since selective angiography of individual collateral vessels, which is tedious and time-consuming, must be performed to ensure success of the treatment;

• non-hepatic arteries that originate from the hepatic artery must be identified to prevent complications such as infarction of other organs.

A catheter is conducted to the branch of the hepatic artery supplying the tumour. To confirm the vascular anatomy and place the catheter tip in the feeder vessel ('superselective catheterization'), rotational Figure 3.17: Three-dimensional imaging for selective chemoembolisation of a liver tumour. (a) Zoomed DSA image. Feeder 1 looks quite easy to catheterize, but it is less clear from which vessel exactly feeder 2 originates. Another concern is the position of the catheter during the chemoembolisation in order to avoid the cystic artery which cannot be identified in the DSA. (b) 3D reconstruction of the anatomy allows for visualization of the cystic artery that is originated from feeder 2. Feeder 2 is branching in a purely anterior direction, which can be seen and understood only from the 3D image. Another key information is that the catheter should be advanced beyond feeder 2 (and not placed in the main vessel) before injection in order not to embolize the cystic artery.

angiography is performed. After which, the treatment can be applied: the chemotherapeutic drug is either emulsified with a iodinated radioopaque oil (Lipiodol by Guerbet) or delivered with calibrated drug-eluting beads as embolic agents, the exact ratio between both substances being adjusted depending on the size and vascularity of the target tumour. The emulsion is concentrated within the liver tumor rather than the surrounding healthy hepatic parenchyma. An interval of several months between each chemoembolisation is necessary for preventing tumour regrowth, while preserving the liver function.

Providing a 3D reconstruction of the liver arterial system removes the anatomical ambiguities and guides superselective catheterization with more accuracy. This is illustrated in the clinical case that is displayed in Fig. 3.17. In particular, the cystic artery often originates off the right hepatic artery and may be missed in a 2D visualization. Arterial angiography is also used for imaging the contrast uptake of the tissues within the liver (liver parenchyma). The interest of rotational angiography for abdominal imaging was demonstrated by clinical studies: a study on 240 hepatic arterial interventions led by Wallace et al. [START_REF] Wallace | Impact of C-arm CT on Hepatic Arterial Interventions for Hepatic Malignancies[END_REF] and a study on 49 patients with hepatocellular carcinoma undergoing chemoembolisation led by Kakeda et al. [START_REF] Kakeda | Usefulness of Cone-Beam Volume CT with Flat Panel Detectors in Conjunction with Catheter Angiography for Transcatheter Arterial Embolization[END_REF] showed that in 41% of the cases and respectively, in 81% of the cases, rotational angiography had impact on patient management. Subtraction of the soft tissues is less important in 3D and is difficult due to respiration. It can be performed under general anaesthesia allowing subtraction of embolised vessels to assess the success of the procedure, that is generally demonstrated by post-procedural CT imaging as illustrated in Fig. 3.19.

Sampling issues in rotational angiography

We have exposed in Sec. 2.1.3 the sampling requirements for tomographic reconstruction. Rotational angiography, however, does not meet these requirements and yields volumes that are deteriorated by sampling artifacts. interventional radiology. It would provide, in particular, straightforward improvements in neuroradiology where collected data are not much affected by motion, but the limited sampling of dense structures (e.g. vessels, coils, metal implants) may generate strong streak artifacts that prevent accurate depiction of the soft-tissues.

Dynamic data issues

We have discussed in Sec. 2.3 the requirements for imaging dynamic data with CT scanners, so that no motion artifacts arise during reconstruction. Because physiological times, e.g. heart beat (about 1 second) or breathing time (range from 3 to 5 seconds) are of the same order of magnitude as the gantry rotation speed of C-arm systems (5 seconds for a 200 • rotation at 40 • per second), cardiac and respiratory motions can affect the performance of rotational angiography. In particular, DSRA is more challenging than DSA, since motion occurring during the protocol (i.e. any of the two scans) would hamper both the reconstruction and the subtraction tasks. The current clinical practice, based on FDK can only ignore motion as the best way to maximize image quality, as was illustrated in Sec. 2.1.3. The relevance of this static assumption however strongly depends on the location of the interventional application:

• in neuroradiology, FDK performs well enough for clinical use [5]: potential displacements of the vessels are very small and localized, and the resulting deterioration is not significant enough to impair the medical interpretation;

• in abdominal imaging, respiratory and cardiac-induced motions frequently arise [START_REF] Van Den Berg | Using three-dimensional rotational angiography for sizing of covered stents[END_REF][START_REF] Tanigawa | Three-dimensional angiography using rotational digital subtraction angiography: usefulness in transarterial embolization of hepatic tumors[END_REF] and image quality may be considerably deteriorated by mild motion blur and strong vessel streak artifacts;

• in cardiac imaging, although not addressed in this work, the temporal variations introduced by the cardiac motion are so important that no static reconstruction is feasible [START_REF] Blondel | Modélisation 3D et 3D+t des artères coronaires à partir de séquences rotationnelles de projections rayons X[END_REF][START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF][START_REF] Neubauer | Clinical feasibility of a fully automated 3D reconstruction of rotational coronary X-ray angiograms[END_REF].

Contrast variations may arise for all the anatomies.

Sampling problems associated with motion

Respiratory motion Respiratory motion affects all the deformable structures in the thorax and abdomen (typically 1 to 2 cycles within a scan of 150 views). Figure 3.21 provides a visual comparison of the FDK reconstruction of the same anatomy in the absence of motion (first row of Fig. 3.21) and in the presence of motion (second row of Fig. 3.21). General anaesthesia ensures patient immobility and allow controlled breath hold, but many interventional procedures are carried out under local anaesthesia only and sedation. To prevent the effect of respiration, the patient is asked to hold his breath during the scanning procedure. Breath hold on a voluntary basis is rendered difficult by sedation and respiratory motion within the scan is often observed in practice. Because only one to two respiratory cycles can be observed within a scan, it is not possible to determine subsampled phases from the acquisition. Our approach casts the reconstruction of dynamic data as a sampling problem that involves the identification of several phases. In this manuscript, we therefore do not deal with the correction of patient motion, including respiratory motion.

Vessel motion Cardiac motion induces small pulsatile displacements of the vascular structures (about 5 cycles within a scan of 150 views), even in organs such as the brain or the liver. When imaging an anatomy that is subject to such motions, it is generally no longer possible to collect uniformly sampled data in time (see Fig. 2.8(b)). Multiple acquisitions synchronized with an ECG would collect uniformly sampled data (see Fig. 2.5(b)), but this is mechanically too demanding (prospective and adaptation of the speed of the gantry to the acquisition).

Contrast variations Besides motion issues, another drawback of lengthy scan times that is specific to angiography is the need for a longer contrast injection to obtain a sufficient and constant enhancement of the blood vessels that are located inside the field of view throughout data acquisition. Contrast variations that arise while scanning yield severe artifacts in the reconstruction. We can distinguish two sources of contrast variations:

• ideally the system is set up with an automated contrast injector that injects the contrast agent at a predefined time point ahead of the scanning. This delay time between the contrast injection and the acquisition is required to fully enhance the blood vessels. Improper synchronization of the injection leads to early or delayed opacification, i.e. the blood vessels appear fully enhanced over a limited angular sector only as illustrated in Fig. 2.8(a);

• in addition, the inadequate injection of the contrast results in incomplete dilution of the contrast dye by the incident blood flow, usually in large vascular structures. The contrast-enhancement varies during the scan. The sampling of fully enhanced blood vessels is cyclic and coarse and non-uniform, similarly to the sampling that is associated to small displacement of the vessels, as shown in Fig. 2.8(b).

Although opacification is a continuous process, we will assume it can be modelled as a binary process of two states 'opacified'/'nonopacified'.

Temporal management

The variations encountered in rotational angiography are so challenging that they cannot usually be corrected for through the motion compensation approaches that were described for CT scanners in Sec. 2.3.4.3. The study of the scan as a temporal sequence is very complex, since the perceived temporal variations within the scan have various causes. To illustrate this discussion, let us analyse a real rotational scan of which a few views were presented in Fig. 3.4. From one 2D view to another, the perceived variations that are due to:

• the motion of the imaging system that is completely described by the projection matrices;

• the magnification of the 3D real object caused by the cone-beam geometry. As explained in 3.1.2 the magnification factor depends on the relative positions of the X-ray source, the imaged object and the detector. Measurements consequently depend on the depth of the imaged structure. Thus, we cannot directly correlate the magnitude of the projected motion in the 2D view to the magnitude of the real 3D motion, since two structures are not affected with the same magnification factor;

• the projection of the three-dimensional motion that affects anatomical structures that are located within the field of view. Because of the projective operation, one dimension of the real motion is lost. However the motion component along the vertical z-axis can be observed in the two-dimensional • the contrast variations in the injected vessels that are located within the field of view.

Let us assume we can decompose the scan into two layers: (a) a first layer that corresponds to the projection of a static background, and (b) a second layer that corresponds to the projection of moving objects (e.g. catheter, organs) with motions of different magnitudes and natures. The perception of the motion is related to the temporal resolution: the higher the frame rate is, the better the motion can be understood, because the temporal resolution is high, and thus the variations between two consecutive views are small. In order to extract the temporal signal that is associated to the vertical component of the catheter motion, we tracked the position along the z-axis of two particular points of the tip (end point and point of inflexion) and plotted this signal in Fig. 3 To explicitly address the temporal variations within the projection data, a classical approach consists in selecting subsets of views from which reconstruction via classical algorithms is feasible. Each time point defines a subset of projections that must ideally fully sample the volume at this particular time point [START_REF] Prümmer | Cardiac C-Arm CT: A Unified Framework for Motion Estimation and Dynamic CT[END_REF]. Since one single rotation is performed only, such an approach presents the main drawback to eliminate a large amount of data and thus to lose an amount of information which is non negligible given the small number of cardiac cycles that are available during a rotational acquisition. Furthermore, synchronization with an ECG is not standard with current C-arm systems and is not applicable to contrast variations, respiratory motion or patient motion of the patient, that are observed with C-arm systems. Chapter 4

Compressed sensing-based iterative reconstruction

We presented in Sec. 2.1.3 a discretization of the inverse Radon transform. As an alternative, iterative strategies translate the tomographic problem into a discrete problem. Then a density mapping is estimated, that fits the projection data by optimizing a function called criterion or cost function (e.g. a quadratic data fidelity term).

Because the obtained solution may not be uniquely determined by the data, the optimized criterion has to be combined with penalties. In particular, the recent developments of the compressed sensing (CS) theory highlight the importance of sparse penalties for handling subsampled data. These penalties are incorporated within a least-square approach through proximal algorithms. Based on the previous works of [START_REF] Capricelli | Algorithmes de projections convexes généralisées et applications en imagerie médicale[END_REF][START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF][START_REF] Pustelnik | Méthodes proximales pour la résolution de problèmes inverses: application à la tomographie par émission de positrons[END_REF], we propose in this chapter a unified framework for the reconstruction of angiographic data. Note that even though weighted least-square terms can be used to model physical degradations and noise, we shall focus our discussion on subsampling only.

Penalized weighted least square (PWLS)

At this stage, it is useful to define the discrete notations and model we will use in the rest of this manuscript. Let us denote f ∈ R K the vector that is associated with the imaged object, where K is the number of voxels in the 3D space. Coefficient (f ) k describes the linear attenuation at voxel k. Let us denote p ∈ R J×N the vector that refers to the full set of measurements, where N is the number of angular positions that are considered by the acquisition and J the number of measurements collected by the 2D detector array at each angle. Coefficient (p) jn represents the log intensity (see Sec. 2.2.2) at pixel j of the detector at angle φ n .

Fidelity term

The relationship between the projection operator and the imaged object is linear. Let us then denote R ∈ R J×N × R K the projection matrix that models the rotational CB acquisition of p from f , where the coefficient (r) k,jn represents the contribution of voxel k to the measure detected at pixel j at angle φ n . Hence, a column of R refers to a given voxel, while a row refers to a given measurement. Let us denote R T the transpose operator of R, in other words, the backprojection operator. The tomographic reconstruction problem consists in calculating an estimation f ∈ R K of f knowing p and R, i.e. solving the system of linear equations:

Rf = p. (4.1)
Note that this formalism is valid for all geometries (parallel-beam, fan-beam and cone-beam) and for all sampling patterns (truncated, uniformly subsampled and nonuniformly subsampled). When projection data are limited in angular range, the system of linear equations of Eq. (4.1) contains fewer equations than unknowns: it is then said underdetermined. On the other hand, it is said overdetermined (or fully sampled) if there are more (independent) equations than unknowns. For exact reconstruction of f , critical sampling requires as many equations as unknowns. To further improve the reconstruction model, it is necessary to incorporate physical degradations (e.g. beam hardening [START_REF] De Man | An iterative maximum-likelihood polychromatic algorithm for CT[END_REF]) or noise [START_REF] Thibault | A Three-dimensional Statistical Approach to Improved Image Quality for Multi-slice Helical CT[END_REF][START_REF] Ramani | A splitting-based iterative algorithm for accelerated statistical X-ray CT reconstruction[END_REF]. Here, however, we will not put forth any model for these nonidealities.

A traditional approach for solving systems of equations such as Eq. (4.1) consists in approximating f by the best fit in the least-square sense, i.e. to find the solution that minimizes the ℓ 2 -norm of the difference between the forward projections and the projection data. It was generalized by Herman [START_REF]On the noise in images produced by computed tomography[END_REF] who introduced a positive-definite matrix W for integration of noise information in the reconstruction model, that is consistent with Bayesian estimation assuming zero-mean Gaussian noise of covariance W -1 : fQ W = arg min

f ∈R K 1 2 Rf -p 2 2,W = Q W (f ) , (4.2)
where Rf -p

2 2,W = (Rf -p) T W (Rf -p).
Provided that vector f is initialized f (0) = 0, the obtained estimate is referred to as the weighted least-square solution of minimum norm that we denote f * and implies inversion of the singular values of R. The explicit solution of Eq. (4.2) is known as the generalized inverse or pseudo-inverse:

f * = (R T W R) -1 R T W p. (4.3)
Note that if f * is such that Rf * = p then this solution does not depend on W . Equation 4.3 requires inverting large size matrices and may be intractable in practice. Because the functional that is minimized in Eq. (4.2) is quadratic, the computation of f * can be handled by gradient descent iterative algorithm:

Gradient descent Initialisation i = 0 f (0) = 0 Algorithm i ≥ 1 f (i) = f (i-1) -τ (i-1) = ∇Q W (f (i-1) ) R T W (Rf (i-1) -p) = (I -τ (i-1) R T W R)
convergence condition

f (i-1) + τ (i-1) R T W p . (4.4)
where I refers to the identity operator and τ (i) > 0 to the gradient step at iteration i.

The estimate at iteration i + 1, f (i+1) , is modified by the addition of the discrete weighted backprojection of the difference between the measurement vector and projection data associated with estimate f (i) , that is scaled by τ (i) . The convergence rate of such an algorithm depends on the norm of operator R T W R and the value of step τ (i) , that is allowed to change at every iteration. The simplest solution consists in choosing constant τ (i) = τ that must be such that I -τ R T W R 2 < 1, but convergence can generally be made faster by using an adaptive gradient step. In the particular case W = I (zero-mean constant variance Gaussian noise) the gradient descent is referred to as the well-known Landweber algorithm. In [START_REF] Payot | Reconstruction vasculaire tridimensionnelle en imagerie par rayons X[END_REF], the convergence rate of the Landweber algorithm has been shown not to be the same for all frequencies, and that in particular low frequencies are recovered faster than high frequencies, hence resulting in a slow convergence rate. A faster alternative to the gradient descent is the conjugate gradient method [START_REF] Culioli | Introduction à l'optimisation[END_REF][START_REF] Kim | An Interior-Point Method for Large-Scale L1-Regularized Least Squares[END_REF] that replaces the gradient direction ∇Q(f (i-1) ) with conjugate descent directions α such that α (i) |R T W Rα (j) = 0 ∀j < i.

Penalty term

We have seen in Sec. 2.1.3 that matrix R T W R is not well-conditioned. The problem is sensitive to noise: for small singular values, noise is large with respect to signal. To improve problem conditioning, apodization can be expressed through (smooth) regularization of Eq. (4.2) [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF]:

fQ W + Af 2 2 λ = arg min f ∈R K Q W (f ) + λ Af 2 2 , (4.5)
where λ sets the relative importance of the regularization with respect to the fidelity term and A is the gradient at a chosen order α, i.e. A = ∇ α . It leads to a trade-off between an unbiased (i.e. a perfect fit to the projection data) but noisy solution and a smooth but biased solution. Note that the impact of the constraint now heavily depends upon the choice of W . The minimization of Eq. (4.5) is carried out using a descent algorithm, similarly to the minimization of Eq. (4.2). To further constraint the reconstruction problem, it may be desirable to introduce penalties that are convex but not differentiable. In that case, however, accelerated gradient descent techniques are no longer applicable.

Proximal algorithms for PWLS

In this section, let us consider the unconstrained minimization problem:

fQ W +χ λ = arg min f ∈R K {Q W (f ) + λ • χ(f )}, (4.6)
where χ is a convex but not necessarily differentiable functional with penalty weight λ. To solve Eq. (4.6), Combettes et al. [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF] propose to split the optimization process so that Q W and χ are individually minimized. The minimization of Q W is thus handled by a gradient descent iteration as described in Sec. 4.1, while the minimization of χ requires introducing new tools.

Proximal Splitting

Let us denote by proj ιE (f ) the projection operator onto convex set E ⊂ R K . It corresponds by definition to the unique minimum of the minimization problem

proj ι E (f ) ≡ arg min g∈R K ι E (g) + 1 2 g -f 2 2 , (4.7)
where ι E refers to the indicator function:

ι E (f ) = 0 if f ∈ E +∞ otherwise .
In particular, let us denote by ι + the indicator function:

ι + (f ) = 0 iff ∈ R K + +∞ otherwise ,
where R K + is defined as the convex set of positive images. When replacing ι E by λ•χ, let us denote by prox λχ (f ) the proximal operator that allows for extending the notion of projection operator:

prox λχ (f ) ≡ arg min g∈R K χ(g) + 1 2λ g -f 2 2 . (4.8)
The forward-backward proximal algorithm builds a sequence of solutions that converges to the minimum of Eq. (4.6) by splitting the minimization into an explicit gradient step for minimizing the data fidelity term and an implicit step applying the constraint through the corresponding proximal operator:

Forward-backward algorithm Initialisation i = 0 f (0) = 0 Algorithm i ≥ 1 f (i) = prox τ λχ implicit step (f (i-1) -τ ∇Q W (f (i-1) ) explicit gradient step
). (4.9)

Fast variants of the forward-backward proximal algorithm, that rely on accelerated gradient methods in the spirit of Nesterov's work in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] are applied to solve problems of the form given in Eq. ( 4.6) at marginal extra computational cost compared to the forward-backward algorithm. Let us briefly present FISTA (Fast Iterative Shrinkage Thresholding Algorithm) that was developed by Beck et al. in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] and that consists in maintaining two variables instead of one at each iteration.

FISTA Initialisation i = 0 f (0) = 0 g (0) = 0 t (0) = 1 Algorithm i ≥ 1              f (i) = prox τ λχ g (i-1) -τ ∇Q W (g (i-1) ) t (i) = 1 + 1 + 4t (i-1) 2 2
g (i) = f (i) + t (i-1) -1 t (i) f (i) -f (i-1) . (4.10)
Such an approach brings the conjugate gradient like super-linear convergence to ℓ 1 -norm penalty. Unlike the simple proximal scheme, there are no convergence guarantees on the sequence of iterates generated by FISTA but has proven to be quite reliable in practice. A monotone variant of FISTA, called MFISTA has been introduced in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

Proximal operator computation

One of the key aspect of proximal splitting is that the proximal operator can be computed either directly or at a reduced computational cost. In the following, we present the computation of some proximal operators that are associated with ℓ 2 -norm and ℓ 1 -norm penalties.

ℓ 2 -norm penalty

Let us consider a ℓ 2 -norm penalty, i.e. χ(f ) = Af 2 2 , where A is a linear operator. The corresponding proximal operator is defined as:

prox λχ (f ) ≡ arg min g∈R K Ag 2 2 + 1 2λ g -f 2 2 = I + λA T A -1 f. (4.11)
This solution can be computed with all gradient techniques for least-squares, that were described in Sec. 4.1.

ℓ 1 -norm penalty

Let us consider a ℓ 1 -norm penalty, i.e. χ(f ) = Af 1 , where A is a linear operator. The corresponding proximal operator is defined as:

prox λχ (f ) ≡ arg min g∈R K Ag 1 + 1 2λ g -f 2 2 . (4.12)
With ℓ 1 -norm, optimization is said to be non-smooth because the gradient of the functional does not exist at the minimum. Depending on the constraint, proximal operator can be computed directly by thresholding or using iterative schemes. We list the key proximal operators we shall use further on, following the presentation of [START_REF] Capricelli | Algorithmes de projections convexes généralisées et applications en imagerie médicale[END_REF][START_REF] Pustelnik | Méthodes proximales pour la résolution de problèmes inverses: application à la tomographie par émission de positrons[END_REF].

Image thresholding

The proximity operator that is associated to the image ℓ 1 -norm (A = I) is the soft-thresholding operator S λ (f ) of threshold λ. Soft-thresholding here can be seen either as a denoising, a compressing or a segmenting step that does not modify the image appearance above the given threshold.

prox λχ (f ) = S λ (f ) =    f + λ if f ≤ -λ 0 if -λ < f < λ f -λ if f ≥ λ (4.13)
Wavelet thresholding Let us denote A w an orthonormal wavelet transform. Among all wavelet transforms, let us cite Haar transform, which is the first and the simplest discrete wavelet transform, or orthogonal Daubechies transform [START_REF] Daubechies | Orthonormal bases of compactly supported wavelet[END_REF]. The proximal operator associated to the wavelet ℓ 1 -norm (A = A w ) is similarly defined as the soft-thresholding operator in the wavelet domain E Aw [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], which can be decomposed into wavelet transform A w , followed by a soft-thresholding operation and inverse wavelet transform

A -1 w : prox χ (f ) = A -1 w [S λ (A w (f ))] . (4.14)
Total variation (TV) penalty To enforce the piecewise constant nature of an image (i.e. to preserve the edges and smooth out other areas), a common approach that was pioneered by Rudin et al. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] relies on minimizing the anisotropic total variation (TV), that is defined for a smooth image as

TV(f ) = ∇f 1 = K-1 k=0 |(∇f ) k |. (4.15)
Chambolle [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF] proposed a dual approach to transform the minimisation of Eq. (4.12) for A = ∇ into a differentiable problem. In such an approach, the ℓ 1 -norm is cast as solution of:

Af 1 = arg max h∈A * (R K ) h ∞ ≤1 Af |h = arg max h∈A * (R K ) h ∞ ≤1 f |A * h (4.16)
where A * denotes the adjoint operator. The TV proximal operator is consequently given by:

prox λχ (f ) = arg max h∈∇ * (R K ) h ∞ ≤1 arg min g∈R K 1 2 g -(f -λ∇ * h) 2 2 - 1 2 f -λ∇ * h 2 2 , (4.17)
where ∇ * corresponds to the opposite of the divergence operator -div. To solve Eq. (4.17), it is possible to use a fixed point iteration method [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF] or a proximal splitting method [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. An alternative implementation of the minimisation of Eq. (4.12) for A = ∇ is the digital TV filter of Chan et al. [START_REF] Chan | The Digital TV Filter and Nonlinear Denoising[END_REF],

that relies on both lagged fixed-point iteration method and regularization |∇u| ǫ = |∇u| 2 + ǫ 2 with ǫ > 0 and consists in linearizing the TV minimization at each iteration.

Nonlinear anisotropic diffusion

The proximal formalism covers all standard image processing operations that are solution of a least-square based variational approach that are mainly used for smoothing and segmentation purposes. In particular, we shall investigate the use of nonlinear anisotropic diffusion [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF] in replacement of the TV penalty:

prox λχ (f ) = arg min g∈R K 1 4 ϕ( ∇g 2 2 ) 2 2 + 1 2λ g -f 2 2 , (4.18)
where ϕ is a real function that locally modulates the diffusion process.

Multiple criterion

So far, we considered that the penalization term captures a single penalty. Let us now consider that it is expressed as a linear combination of penalties.

Indicator functions Let us first consider the particular case of a multiple criterion that is a combination of indicator functions:

χ(f ) = C c=1 ι Ec (f ), (4.19)
where E c ⊂ R K . Minimization of Eq. (4. [START_REF] Candès | Decoding by linear programming[END_REF]) is carried out through successive activation of each convex set:

f (i) = C c=1 proj ι Ec (f (i-1) ). (4.20)
This process is known as projection onto convex sets (POCS) [START_REF] Youla | Image Restoration by the Method of Convex Projections: part 1. Theory[END_REF]. Many reconstruction algorithms derive from POCS, among which the well-known algebraic reconstruction technique (ART) [START_REF] Gordon | Algebraic Reconstruction Techniques (ART) for threedimensional electron microscopy and x-ray photography[END_REF] that converts the least-square data fitting term into a set of penalties, one per measurement.

Image thresholding with positivity Let us combine the image ℓ 1 -norm with the projection on the convex set of positive images, which is separable with respect to the image domain:

χ(f ) = λ • f 1 + ι + (f ). (4.21)
Because the image ℓ 1 -norm is also separable, it has been shown in [START_REF] Pustelnik | Méthodes proximales pour la résolution de problèmes inverses: application à la tomographie par émission de positrons[END_REF] that the resulting proximal operator is simply defined by:

prox χ (f ) = proj ι+ •S λ (f ) = f -λ if f ≥ λ 0 otherwise . (4.22)
Note, however, that this direct computation of the proximal operator does not hold for arbitrary penalty χ and arbitrary convex set E ⊂ R K .

Generic case

Let us now consider the following multiple criterion:

λ • χ(f ) = C c=1 λ c • χ c (f ) (4.23)
The minimization of the proximal operator associated to Eq. (4.23) can be achieved by the parallel Dykstra-like proximal algorithm [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF], that is given in pseudocode hereafter, provided that prox λcχc (f ) is available and

C c=1 λ c = 1.
Parallel Dykstra-like proximal algorithm

Initialisation i ′ = 0 v (0) = u g (0) 1 = v (0) . . . g (0) C = v (0) Algorithm i ′ ≥ 1                                                 h (i ′ ) 1 = prox τ χ1 g (i ′ -1) 1 . . . h (i ′ ) C = prox τ χ C g (i ′ -1) C v (i ′ ) = C c=1 λ c • h (i ′ ) c        g (i ′ ) 1 = v (i ′ ) + g (i ′ -1) 1 -h (i ′ ) 1 . . . g (i ′ ) C = v (i ′ ) + g (i ′ -1) C -h (i ′ ) C (4.24)

Compressed sensing (CS)

Let us now shortly recall the sampling theory discussed in Sec. 2.1.3: the Shannon-Nyquist criterion provides a sufficient sampling condition for accurate reconstruction, but not a necessary one. If further restrictions are imposed on the signal, then less stringent conditions may be sufficient.

Sparse optimization

The recent developments of the mathematical theory of CS provides a sub-Nyquist sampling criterion that basically depends on signal sparsity: intuitively, the sparser a signal is, the fewer data samples should be required to obtain an accurate reconstruction.

ℓ 0 -pseudo norm penalty The 'best' measure for sparsity is arguably defined as the ℓ 0 -pseudo norm of f which corresponds to the cardinality of the support of f in some transform domain. A signal f is said S-sparse if there exists a domain in which the signal f of size K can be expressed with S non-zero coefficients only (such that S ≪ K):

S = Card{(Af ) k = 0| 1≤k≤K } ≡ Af 0 (4.25)
where Card{•} refers to the cardinal of a set and • 0 denotes the ℓ 0 -pseudo norm. Recovery of f is obtained by identifying the sparsest solution among the solutions of the system of equations presented in Eq. (4.1) through resolution of the combinatorial optimization problem:

fCS-ℓ0 = argmin

f ∈R K Af 0 s.t. Rf = p. (4.26)
Unfortunately, ℓ 0 -minimization is not a convex optimization problem, and the problem described in Eq. (4.26) is computationally intractable, except for some particular cases.

Sparse norm penalties

The ℓ 0 -pseudo norm is fortunately not the only norm that promotes sparsity. The ability of a norm to promote sparsity is related to the geometry of the unit ball which is associated to the norm. This is illustrated in Fig. 4.1 that displays a space with two degrees of freedom, in which the collection of measurements that satisfy Eq. (4.1) is described as a line, and that compares the minimum ℓ 2 -norm solution to the minimum ℓ 1 -norm solution. The unit ball that is associated to the ℓ 2 -norm is isotropic and does not favor any specific direction of the space as shown in Fig. 4.1(a). The minimum ℓ 2 -norm solution is obtained by taking the circle around the origin that is tangent to the solution line. Large coefficients are therefore penalized more heavily and the solution is generally not sparse. On the other hand, the unit ball that is associated to the ℓ 1 -norm is anisotropic: it is a square with singular points due to the discontinuities of the ℓ 1 -constraint, that are located along the axis of the coordinate system, as shown in Fig. 4.1(b). The minimum ℓ 1 -norm solution happens to intersect the line at one of those points. Because, it has only one coordinate different from zero, it is sparser than the minimum ℓ 2 -norm solution. Many small coefficients tend to carry a larger penalty than a few large coefficients and are consequently suppressed. Note that when focusing on the cardinality of the solution only, no side information about the patterns of relevant (non-zero) coefficients is provided. As an alternative, it is possible to consider a structured sparsity to simultaneously select all variables forming a group. Such an approach requires penalizing other functions than the cardinality of the support or regularizing by other norms as suggested in [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF]. Candès theorem As shown by Candès et al. in [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF], CS states that if a signal can be made sparse or compressible in some transform domain that is incoherent with the space in which the sampling is performed, then it is possible to reconstruct the signal accurately from fewer samples. Exact recovery of f is obtained by solving Eq. (4.26). In [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF], a simple way is proposed to make the problem described in Eq. (4.26) tractable, that consists in casting the sparse estimation as a convex optimization problem by replacing the ℓ 0 -pseudo norm by a ℓ 1 -norm:

0 p = Rf f 2 = 1 0 p = Rf f 1 = 1
fCS-ℓ1 = argmin f ∈R K Af 1 s.t. Rf = p. (4.27)
The key result of the CS theory is that the ℓ 1 -minimizer of Eq. (4.27) is unique, and leads to the sparsest solution, i.e. solutions of Eq. (4.26) and Eq. (4.27) are equivalent. An upper bound S max for the maximum level of sparsity for which exact reconstruction can be achieved is proposed in [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF] under the following assumptions:

• a sparse representation of the solution in a known transform domain exists;

• the sampling is random;

• matrix R obeys the restricted isometry property (RIP) for δ 2S < √ 2 -1 [START_REF] Candès | Decoding by linear programming[END_REF] where δ S denotes the isometry constant of R, that is defined as the smallest number such that for all S-sparse vectors:

(1 -δ S ) • f 2 2 ≤ Rf 2 2 ≤ (1 + δ S ) • f 2 2 .
Interestingly, Candès et al. [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF] illustrate the CS theory by showing that the exact reconstruction of the Shepp-Logan phantom test image -which is a piecewise constant image-is achievable from a few amount of integral lines only through minimization of the image TV even though the random sampling and RIP property are not verified for R. This result opened tremendous perspectives for image quality improvement in tomography.

Greedy algorithms for ℓ 0 -pseudo norm minimization

Some works propose a heuristic approach known as Orthogonal Matching Pursuit (OMP) to build the CSsolution fCS . Orthogonal matching pursuit belongs to the iterative greedy methods that build a sequence of locally optimal optimizations to determine a globally optimal solution. The major advantage of this type of approach is that it admits simple, fast implementations [START_REF] Mallat | Matching Pursuit with time-frequency dictionaries[END_REF][START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Mallat | Adaptive time-frequency decompositions[END_REF][START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF].

An intuitive illustration that OMP recovers the image support in tomography is provided in Fig. 4.2. It is based on thresholding the structure having the highest intensities in the reconstructed image and then apply R to these structures to subtract their contribution to the real measurements. The thresholding also removes the subsampling pattern that is associated to the thresholded structures, and hence recovers lower intensity structures that were previously hidden by the artifacts. By iteratively repeating this procedure, it is possible to recover all the structures of a sparse image. This approach was proposed for MR images in [START_REF] Fain | Correction for artifacts in 3D angularly undersampled MR projection reconstruction[END_REF]. A similar procedure was proposed by Donoho et al. [START_REF] Donoho | Sparse Solution of Underdetermined Systems of Linear Equations by Stagewise Orthogonal Matching Pursuit[END_REF], called Stagewise Orthogonal Matching Pursuit (StOMP). Let us index by support Ω in the image domain, the restricted application of the (back-)projection operator to Ω. In StOMP, a number of iterations or stages is fixed a priori. At each stage s, an estimate of the reconstruction support is selected, that is best correlated with the residual Rf (s-1) -p, this support is then merged with the previously estimated support Ω s-1 to form the updated support Ω s . StOMP finally performs least-square reconstruction on the limited optimal support Ω = Ω s :

fStOMP = arg min f ∈Ω 1 2 R Ω f -p 2 2,W . (4.28)
Note that when the data are not sparse, Ω tends to the entire space, and hence this approach results in the standard least-square reconstruction.

Basis pursuit

The minimization of Eq. (4.27) is known in the literature as basis pursuit (BP) [START_REF] Chen | Atomic Decomposition via Basis Pursuit[END_REF] and remains a preferred method to promote sparsity even though ℓ 1 penalized solutions are not always sparse, they need not be unique in general, and they need not lead to the sparsest solution. One of the main advantages of BP is that it enjoys rigorous proofs of exact reconstruction under quite general circumstances. Note that penalizations by ℓ p -norms with p < 1 have also been proposed, but are non-convex problems, for which finding global minimizers is more tedious. Optimization methods adapted for solving the BP problem described by Eq. (4.27) have been a focus of research interest before CS was first introduced. Traditional optimization tools rely on a necessary and sufficient optimality criterion that is valid for differentiable functionals only. Because the ℓ 1 -norm is not differentiable, the approaches that are developed in the literature for solving Eq. (4.27) include:

• regularization of the ℓ 1 -norm so that it becomes differentiable [START_REF] Chan | The Digital TV Filter and Nonlinear Denoising[END_REF] but it breaks to some extent the capacity of the penalty to induce sparsity; • casting of the BP problem in a standard linear program which is convex and tractable. It is then solved using general purposed solvers such as simplex and interior point methods [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Kim | An Interior-Point Method for Large-Scale L1-Regularized Least Squares[END_REF] for small scale problems and proximal algorithms for large scale problems. An alternative linearisation is proposed in [START_REF] Figueiredo | Gradient Projection for Sparse Reconstruction: Application to Compressed Sensing and Other Inverse Problems[END_REF];

• subgradient descent that can be used as soon as a subgradient can be computed efficiently [START_REF] Shor | Minimization methods for non-differentiable functions[END_REF]. This last strategy is widely applicable, has low running time complexity per iterations, but has low convergence rate and ususally yields non-sparse solutions.

These criteria are also robust with respect to noise: basis pursuit (BP) will still recover a good approximation f to f if Eq. (4.27) is casted as an ℓ 1 -minimization with a quadratic inequality constraint:

f * = argmin f ∈R K Af 1 s.t. Q W (f ) ≤ ǫ, (4.29)
where ǫ is a nonnegative parameters. The minimization carried out by Eq. (4.29) is closely related to the Least Absolute Shrinkage and Selection Operator (LASSO) that was extensively studied by the Statistics community [START_REF] Tibshirani | Regression Shrinkage and Selection via the LASSO[END_REF][START_REF] Efron | Least angle regression[END_REF] and that consists in a quadratic optimization with a ℓ 1 inequality constraint:

f * = argmin f ∈R K Q W (f ) s.t. Af 1 ≤ η (4.30)
where η is another nonnegative parameters. Note that there exist in general no direct correspondence between parameters ǫ and η, except for some special cases. Therefore algorithms that are specific for solving one of the two formulations in Eq. (4.29) and Eq. (4.30) may not be used directly for solving the other.

There exists an equivalent non-constrained formulation of Eq. (4.29) and Eq. (4.30) that is a ℓ 1regularized least square:

f * = argmin f ∈R K {Q W (f ) + λ Af 1 } (4.31)
where λ governs the degree of sparsity of the minimizer. All formulations easily integrate multiple sparsity constraints into the reconstruction through a penalization that is expressed as a linear combination of C regularization terms as described in Eq. (4.23). It can be solved using a nonlinear conjugate gradient descent algorithm with backtracking line search [START_REF]Sparse MRI[END_REF]. In this manuscript, we rather rely on proximal splitting methods of Sec. 4.2.

f * = argmin f ∈R K Q W (f ) + C c=1 λ c A c f 1 (4.32)

Homotopy for ℓ 1

Although shown to perform more efficiently than BP in certain circumstances [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF], OMP fails to find the sparsest solution in others where ℓ 1 -minimization succeeds and does not offer in general the strong theoretical guarantees that are attached to ℓ 1 -minimization. In [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF], Tropp provides a sufficient condition under which both OMP and BP can recover the optimal representation of an exactly sparse signal. Moreover, he proves that OMP may be able to recover the support from the optimal representation even when the signal is not perfectly sparse, which is the case of natural signals. A parallel can be established between OMP approaches and homotopy approaches, of which some bridging elements are discussed in [START_REF] Donoho | Fast Solution of ℓ 1 -norm Minimization Problems When the Solution May be Sparse[END_REF]. The key idea behind homotopy is that when varying the regularization hyperparameter λ over (0, +∞), the family of solutions ( fℓ1 ) λ of Eq. (4.31) is defined as the regularization path. Efron et al. [START_REF] Efron | Least angle regression[END_REF] show that there exists a series of values λ 1 , • • • , λ s , with 0 = λ S < • • • < λ 1 = λ max , such that the regularization path is a piecewise linear function of λ:

( fℓ1 ) λ = λ s -λ λ s -λ s+1 ( fℓ1 ) λs+1 + λ -λ s+1 λ s -λ s+1 ( fℓ1 ) λs , λ s+1 ≥ λ ≥ λ s , s = 1, • • • , S,
where ( fℓ1 ) λs is solution of Eq. (4.31) when λ = λ s . To compute the regularization path, homotopy methods were developed by varying the regularization parameter λ from large to small values [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Efron | Least angle regression[END_REF] and proved to be as efficient as greedy stepwise algorithms such as OMP when the solution is very sparse. The number of stages S is a key parameter. It is shown to be equal to the sparsity of the target solution for some special cases in [START_REF] Donoho | Fast Solution of ℓ 1 -norm Minimization Problems When the Solution May be Sparse[END_REF]. If the intermediate solutions are cheap to compute, then they provide a sequence of convenient estimations for initialization to the next problem. In practice, they are only approximated, making them cheaper to compute. Provided that an approximate solution is obtained at the end of each stage, the recent work of Xiao et al. [START_REF] Xiao | A Proximal-Gradient Homotopy Method for the Sparse Least-Squares Problem[END_REF] demonstrates that for large-scale problems, the number of stages can be much smaller than the sparsity the target solution and ensures convergence of the process provided that each problem is solved with sufficient precision and hyperparameter λ is geometrically decreased.

Reconstruction for subsampled data in rotational angiography

Compressed sensing was applied to tomographic, divergent beam problems for the first time by Sidky et al. in [START_REF] Sidky | Accurate Image Reconstruction from Few-Views and Limited-Angle Data in Divergent-Beam CT[END_REF][START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF] that proposed for accurately solving constrained TV-minimization. As discussed by Theriault et al. in [START_REF] Theriault Lauzier | Prior image constrained compressed sensing: Implementation and performance evaluation[END_REF], we must emphasize that straightforward application of the mathematical conditions introduced by the CS framework is not possible in rotational angiography, since it is very challenging, if not impossible, to design a data acquisition method that completely satisfies the mathematical conditions required by the CS theory (see Sec. 4.3).

In this manuscript, we focus on the empirical application of the CS method without pursuing the mathematical rigor of the CS theory itself. In particular, we relax the assumption that is used to establish CS upper bound S max on whether the data are randomly sampled. This is a key difference when applying CS to X-ray imaging with respect to MRI [START_REF] Lustig | Compressed sensing MRI[END_REF][START_REF] Vaswani | Kalman filtered compressed sensing[END_REF]. On one hand, sampling in MRI imaging is a special case of the CS theory: the sampled linear combinations are simply Fourier coefficients and to some extent, a subset of coefficients can be selected. On the other hand, the angular sampling of the spectral domain is imposed by the X-ray CT system. In the context of subsampled X-ray data, sampling issues have been shown to be dealt with through original approaches. Motivated by the desire to reduce dose, Sidky et al. proposed an iterative reconstruction for tomosynthesis [START_REF] Sidky | Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved imagereconstruction algorithms[END_REF] and high resolution image reconstruction [START_REF] Sidky | High resolution image reconstruction with constrained, total-variation minimization[END_REF]. Based on this approach, many works propose to filter out streaks assuming that the object image is sparse in a certain domain.

Note that the validity of compressed sensing is established with least square data fit only. State-ofthe art CS reconstructions mostly rely on an adaptive steepest variant of POCS known as ASD-POCS [START_REF] Sidky | Accurate Image Reconstruction from Few-Views and Limited-Angle Data in Divergent-Beam CT[END_REF][START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF] or ART [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF]. The proximal splitting we presented in Sec. 4.2 actually covers a wide range of published reconstruction algorithms as well as all the algorithms used in this work. Sparse assumptions are generally not valid for clinical images, so that only sparse approximation of the reconstruction problem solution can be obtained. Object sparsity in the image domain was mainly promoted for angiography [START_REF] Li | Improved iterative algorithm for sparse object reconstruction and its performance evaluation with micro-CT data[END_REF][START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF]. Most works focused on piecewise constant assumption because it relies on total variation (TV) [START_REF] Persson | Total variation norm for three-dimensional iterative reconstruction in limited view angle tomography[END_REF][START_REF] Song | Sparseness prior based iterative image reconstruction for retrospectively gated cardiac micro-CT[END_REF][START_REF] Bergner | A comparison of 4D cone-beam CT algorithms for slowly rotating scanners[END_REF][START_REF] Choi | Compressed sensing based cone-beam computed tomography reconstruction with a first-order method[END_REF][START_REF] Ritschl | Improved Total Variation-Based CT Image Reconstruction Applied to Clinical Data[END_REF][START_REF] Sidky | A constrained, total-variation minimization algorithm for low-intensity x-ray CT[END_REF], for which efficient minimization algorithms have been proposed [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Medical images are only approximately piecewise constant, which yields 'patchy' artifacts and unnatural image appearance. Sparsity of the image wavelet transform may be clinically more acceptable, but it has not been proven to be efficient enough alone and is then often used jointly with TV [START_REF] Lustig | Compressed sensing MRI[END_REF][START_REF]Sparse MRI[END_REF].

Compressed sensing with iterative FBP

In this work, we focus on iterative FBP (iFBP), that is referred to as iFDK for the particular case of CB geometry.

Iterative FBP

Mathematically, iterative FBP consists in minimizing functional Q W defined in Eq. (4.2) where matrix W is replaced by the ramp filter D that was introduced in Sec. 2.1.2.2 and is positive and diagonal in the Fourier domain: fiFBP = arg min

f ∈R K Rf -p 2 2,D (4.33)
The iterative scheme of iFBP is described by Eq. (4.4) with W = D and keeping the gradient step τ constant over all iterations, i.e. τ (i) = τ . The solution at iteration i, f (i) , is modified by the addition of the discrete FBP of the difference between the measurement vector and projection data associated with the present estimate, that is scaled by τ .

When the problem is well-determined, matrix R T DR is close to the identity operator and a constant step τ ≈ 1 provides fast convergence. The kernel of the projection matrix, denoted Ker(R), is reduced to singleton {0} and f = f * is the unique solution. For underdetermined problem with truncation, it has been used in conjunction to an ℓ 2 -norm regularization to improve the reconstruction of truncated area. Iterative FDK has also been proposed for correcting for CB artifacts but it has not been shown to converge quickly in that case.

This work however aims at solving the reconstruction problem in the underdetermined case where the measurements are severely subsampled and consequently, the dimension of p is small with respect to the dimension of f . In this case, matrix R T DR is not equal to the identity. Gradient methods may require a smaller step τ ≪ 1. The dimension of Ker(R) is far greater than 1 and f is restricted to an affine subspace of R K , that does not enable complete determination of f . Hence, there exists an infinity of solutions that are compatible with Eq. ( 4.1) and can be decomposed as follows:

f = f * + f ⊥ , (4.34)
where f ⊥ ∈ Ker(R) refers to the component of f that is unobservable with respect to the acquisition sampling. By construction, the unobservable component of initialization f (0) is preserved. In particular, f * sets to zero the unobserved spectral coefficients, which results in subsampling streak artifacts similar to the ones described in Sec. 2.1.3 for the FBP reconstruction. From there it is clear that all images that satisfy Eq. (4.34) are not necessarily relevant from a physical point of view, with impact on image quality.

To reduce the magnitude of f ⊥ , which can be interpreted as the error image made of the complementary streak artifacts, the selection of a clinically relevant solutions is carried out through the introduction of constraints. This criterion can be minimized with proximal algorithms. The imaged densities are naturally positive values. Consequently, a positivity constraint can be introduced at a negligible computational cost. The minimization process identifies a positive image that best fits the projection data [START_REF] Langet | Compressed Sensing Based 3D Tomographic Reconstruction for Rotational Angiography[END_REF]:

fiFBP+ = arg min f ∈R K {Q D (f ) + ι + (f )}. (4.35)
The FBP and iFBP reconstructions yield similar results only if the tomographic problem is well determined -in other words, if the sampling rate is high. In the underdetermined case, FBP reconstruction result does not satisfy Eq. (4.34). To illustrate this point, we compared the FDK and iFDK reconstructions of clinical data with different subsampling factors. Figure 4.3 shows the reconstruction results on an axial slice that presents air, soft-tissues and bones. To enhance the display of the angular subsampling artifacts, we set WL=1000 HU and WW=2000 HU. For the sampling rate of current C-arm systems (minimum 150 views), FDK and iFDK that are compared in Fig. 4.3(a)) appear very similar, iFDK providing however a slightly higher resolution. When the number of views is reduced as shown in Fig. 4.3(b) for a subsampling factor of 2 to (d) for a subsampling factor of 8, standard FDK contains high-frequency -and hence streak artifacts-due to the ramp. On the other hand, iFDK reconstruction better preserved the information at coarse scales, but lost most high-frequencies such as edges. This difference between FDK and iFDK is significant at very low sampling rates (Fig. 4.3(d)). Some intensity profiles through the axial slice of Fig. 4.3 are plotted in Fig. 4.4. We observe that the subsampled iFDK profile is closer than the subsampled FDK profile (Fig. 4.4(b)) to the two profiles obtained with 150 views (Fig. 4.4(a)). This demonstrates a higher robustness of the iFDK reconstruction with respect to the FDK reconstruction.

Iterative FBP with an image ℓ 1 -norm penalty

Considering the contrast-enhanced vessels that cover only a small fraction of the subtracted volume, the notion of sparsity is straightforward, as illustrated by the MIP representation in Fig. 4.5(1.a). The few amount of non-zero voxels (61% of the voxels have values below 100 HU) is confirmed by the study of the intensity histogram which is displayed Fig. 4.5(1.b). The first clinical application of rotational angiography with C-arm systems relied on this feature to propose a reconstruction of the subtracted 3D vessels for neuro-interventional radiology with no reconstruction of the background structure that were removed by digital subtraction within the projections. Non-subtracted volumes and even subtracted volume containing not only the contrast-enhanced vessels, but also perfused structures such as parenchyma cannot, in general, be assumed sparse in the image domain as it is formulated in Eq. (4.25). Thus, to address the reconstruction of the coronary arteries, Hansis et al. [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF] suppress the background structures via a morphological top-hat filter prior to the reconstruction. A closer look at the intensity histogram calculated from reconstructed angiographic data and displayed in Fig. 4.5.

(2) demonstrates a hierarchical sparsity with respect to the intensity level (1.4% of the voxels with value above 2000 HU while about 76% with value between 900 and 1600 HU) due to the selective contrast-injection into the vessels. In the following, we describe this feature of angiographic data as 'sparse vessels over a non-sparse background'.

As discussed in Sec. 2.1.3, subsampling introduces a pattern of positive and negative streaks in the reconstruction. When dealing with subtracted data, 'dark' streaks mostly take negative values and hence, it is possible to enforce the positivity of the solution in order to mitigate the global artifact pattern. Because the background level is generally not equal to zero, negative streaks result in abnormally lower, yet not necessarily negative, pixel intensities. Although still relevant from a physical point of view, the positivity constraint, no longer provides effective artifact reduction.

Combining the image ℓ 1 -norm with the positivity constraint (cf. Sec. 4.2) promotes sparsity by removing at each iteration all structures whose intensity is lower than threshold λτ , a process that we refer to as soft background subtraction in the following of this manuscript. From the image processing point of view, it can be interpreted as an intensity-based segmentation task. Considering the case of a sparse image that only contains vessels, there is no strictly positive λ value that may not remove some part of the vessels, and thus biase the reconstruction result. Consequently, CS-sbs approach proposed in [START_REF] Langet | Compressed Sensing Based 3D Tomographic Reconstruction for Rotational Angiography[END_REF] defines a set of S decreasing λ sbs values:

Λ sbs = λ (s) sbs |s = 1, • • • , S such that λ (1) sbs ≥ • • • ≥ λ (S) sbs = 0,
and solving the sequence defined by the S corresponding ℓ 1 -regularized problems that are of the form: sbs identifies a sparse approximation in the image space that best fits the data with a level of sparsity that is proportional to λ (s) sbs . In practice, the iteration described by Eq. (4.36) results in a standard FBP reconstruction that is scaled by τ and segmented with threshold τ λ (s) sbs . Thus, CS -sbs is compatible with the streak-free reconstruction of sparse structures first (e.g. high-intensity vessels), while the non-sparse background is progressively reintroduced as λ (s) sbs tends to zero. We observed convergence of the global algorithm when solving each problem described by Eq. (4.36) with one iteration and using a simple linear decrease of λ (s) sbs providing:

fχ sbs λ sbs = arg min f ∈R K Q D (f ) + λ sbs f 1 + ι + (f ) = λ sbs • χ sbs (f ) , ( 4
• λ (1)
sbs is initialized equal to 90 % of the maximum value of the FDK reconstruction,

• while λ (S) sbs ≥ 0 is chosen as small as desired.

Iterative FBP with TV penalty

Assuming the most relevant information in the images is provided by the edges, the image can be sparsified by applying a discrete gradient operation. This is illustrated in Fig. 4.6 that displays the gradient image of a slice and the corresponding intensity histogram.

Minimizing TV has been extensively studied in tomographic reconstruction for the correction of subsampling artifacts [START_REF] Ritschl | Improved Total Variation-Based CT Image Reconstruction Applied to Clinical Data[END_REF][START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF][START_REF] Persson | Total variation norm for three-dimensional iterative reconstruction in limited view angle tomography[END_REF]. These methods show streak artifact reduction, but at the expense of an overall change in the image appearance, that may not be clinically acceptable. Considering a coarse-to-fine approach, we propose in CS-tv to solve a sequence of S ℓ 1 -regularized problems where the constraint is the combination of a constraint on the image gradient ℓ 1 -norm ∇f 1 with the positivity constraint on the image domain: piecewise-constant approximation that best fits the data with a level of sparsity that is proportional to λ (s)

fχtv λtv = arg min f ∈R K Q D (f ) + λ tv ∇f 1 + ι + (f ) = λ tv • χ tv (f ) , ( 4
tv . Total variation minimization as described in Sec. 4.2.2.2. Because computing the gradient of an image is a high-pass filtering operation, the gradient transform can be considered as computing some sort of fine-scale wavelet transform (without computing coarser scales). Performing CS -tv is compatible with the reconstruction of rough piecewise constant approximation first, since strong λ (s) tv value forces the reconstruction result to have smaller total variation (i.e. only very intense structures or large background area are kept), while the fine-scale information is progressively reintroduced as λ (s) tv is decreased. We observed convergence of the global algorithm when solving each problem with one iteration and using a simple linear decrease of λ (s) tv providing:

• λ (1)
tv is initialized so that only coarse structures are kept from FDK reconstruction,

• while λ (S)
tv is either set equal to 0 (i.e. the final solution is unbiased) or tuned so that most of the fine details remain (i.e. the final solution is slightly regularized).

Iterative FBP with wavelet ℓ 1 -norm penalty

The use of wavelet was previously investigated for tomographic reconstruction in MRI by Lustig et al. in [START_REF] Lustig | Compressed sensing MRI[END_REF][START_REF]Sparse MRI[END_REF]. The exploitation of such sparsifying transforms was motivated by the success of image compression that relies on such transforms to reduce data redundancy by encoding the few significant coefficients and storing them, for later decoding and reconstruction with little or no visual loss of information. Each wavelet coefficient carries both spatial position and spatial frequency information. In particular, coarse-scale wavelet coefficients determine the low resolution image components, while fine-scale wavelet coefficients determine the high resolution components. Note that the recovery of coarse-scale wavelet coefficients requires less angles (see Sec. 2.1.3). Koff et al. have shown in [START_REF] Koff | An Overview of Digital Compression of Medical Images: Can We Use Lossy Compression in Radiology?[END_REF] that, to some extent, medical images are compressible, and hence a wavelet transform can be applied without loss of crucial information.

Image compression using Daubechies 9/7 wavelet transform, that is used in JPEG 2000 -one of the most well-known image compression standards-, is illustrated in Fig. 4.7: no visual difference is perceptible between the original reference image and its compression through removal of 90% of the coefficients of the first detail subband, while only slight deteriorations appear for its compression through removal of 90% of the coefficients of the first two detail subbands. Considering a coarse-to-fine approach, we pro-pose in CS-w to solve a sequence of S ℓ 1 -regularized problems where the constraint is the combination of a constraint on the image wavelet transform ℓ 1 -norm Ψ w f 1 , where Ψ w denotes a discrete wavelet transform, with the positivity constraint on the image domain: The image wavelet transform Ψ w can be for instance Haar transform or orthogonal Daubechies transform. As given in Sec. 4.2, the proximal operator is soft-thresholding in the wavelet domain [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF]. The minimization at λ w = λ (s) w identifies a compressible (in the sense of the wavelet transform) approximation that best fits the data with a level of sparsity in the wavelet domain that is proportional to λ (s) w . Note that in medical imaging, standard decomposition is limited to a few subbands only, since compression above three levels of fine-scale subbands raises both legal and diagnostic issues. We observed convergence of the global algorithm when solving each problem with one iteration and using a simple linear decrease of λ (s) w providing:

fχw λw = arg min f ∈R K Q D (f ) + λ w Ψ w f 1 + ι + (f ) = λ w • χ w (f ) , ( 4
• λ (1)
w is initialized so that 90 % of the coefficients in the detail subbands are set to zero in the wavelet transform of FDK reconstruction,

• while λ (S)
w is either set equal to 0 (i.e. the final solution is unbiased) or tuned so that 10 % of the coefficients in the detail subbands are set to zero in the wavelet transform of FDK reconstruction (i.e. the final solution is regularized).

Numerical simulations with uniform subsampling

To establish the behaviour of the CS reconstruction algorithms and, in particular, to evaluate the image quality that can be achieved with respect to a ground-truth, we built a numerical phantom with simulated contrast-enhanced vessels by adding a disk pattern (with intensities from 1500 to 6000 HU) to a 512 × 512 head phantom cross-section where background structures are valued between 1000 and 2000 HU. . We simulated the rotational acquisition of projection data in parallel geometry by integrating the density of the phantom along parallel lines. Note that the effect of photon statistics was not simulated but could be added for more realism.

Homotopy vs. regularization

We simulated the acquisition of projection data with 150 equiangularly distributed views over 180 • . We carried out regularized reconstructions with TV minimization and thresholding in the wavelet domain.

To generate these reconstructions, we used penalized iFBP for hyperparameters λ w and λ tv fixed values.

In addition, we performed reconstruction with homotopy with soft background subtraction, TV minimization and thresholding in the wavelet domain. Reconstruction settings are indicated in Tab. 4.1 and Tab. 4.2.

Parameter Value τ 0.9 S 25 backgrounds I 1 iteration / background Table 4.1: Generic parameters common to all iFBP-derived algorithms.

Figure 4.9(a) shows the reconstruction result for iFBP+. It serves as reference to which compare regularized and homotopy reconstructions. Although minor streaks that originated from the high-intense structures (e.g. simulated contrast-enhanced vessels and skull bone) are visible near the periphery of the skull as highlighted by the detail images, the entire image is relatively free of streak artifacts.

Regularization There are interesting observations that one can make regarding the regularized reconstructions. The reconstruction results shown in Fig. 4.9(b) and in Fig. 4.9(d) correspond to the minimization of a functional that enforces sparsity, TV for experiment n • 2 and respectively wavelet ℓ 1 -norm for experiment n • 4. Sparse penalization removed the thin streak artifacts but at the expense of a significant change in the image appearance: total variation resulted in a piece-wise constant image with a cartoon look that may not be acceptable for clinical purpose, while thresholding in the wavelet domain yielded local artifacts. Such results do not fit the projection data, but are instead the regularized approximation with a level of sparsity that is given by the regularization hyperparameter that best fits the data. We then attributed a weaker weight to the sparse penalization. This corresponds to the reconstruction results with TV penalty in experiment n • 3, that is shown in Fig. 4.9(c), and with wavelet ℓ 1 -norm penalty in experiment n • 5, that is shown in Fig. 4.9(e). These reconstructions appear very close to the iFBP+ reconstruction, and hence the sparse penalization did not have a significant impact on image quality. Consequently, there is a trade-off between sparse regularization and fidelity to the data. Finding the best trade-off is not a trivial task, in particular because it generally depends on the imaged structures. Homotopy The previous study of the iFBP+ reconstruction shown in Fig. 4.9(a) made clear that we are trying to extract more information from relatively finely sampled data. Consequently, the impact of CS is expected to be small. Nevertheless, the use of different sparsity constraints with homotopy, of which the resulting images are compared in Fig. 4.10, provided some visible advantages over iFBP+ and the regularized reconstructions. CS-sbs yielded a non-biased reconstruction that is shown in Fig. 4.10a.

A careful study of the top detail image indicates an orientation change of the streaks that are visible below the skull and that originated from the contrast-enhanced vessels in Figure 4.9(a). The latter were corrected by applying soft background subtraction. The subsampling streaks we now observe were more likely caused by the skull, given their orientation. We used CS-tv and CS-w to generate non-biased (experiments n reconstructions, since mitigation of the streak artifacts is achieved without introducing unrealistic patches or wavelet decoding artifacts. These observations are confirmed by lower RMSD values, provided in Table 4.2. Note that non-biased CS-w did not correct for the subsampling streaks which were corrected by CSsbs, but rather for high frequency streaks, regardless of their origin. Visual inspection of the numerical simulations would lead us to favor CS-tv that targets a slight regularized fit to the data, rather than an unbiased fit. In addition, note that the importance of a obtaining a regularized solution should be reinforced by the presence of noise in the clinical data, a parameter that we did not simulate here.

Homotopy reconstruction for a low uniform sampling rate

We simulated the acquisition of projection data delivering 30 equiangularly distributed views over 180 • . Existing regularization strategies result in a trade-off between a fitting term and a penalty. With 150 views, in the presence of very dense structures (e.g. metal implants, contrast-enhanced vessels), we showed that the CS methodology allows for a more important mitigation of the streak artifacts than regularization, while preserving FBP-like appearance. When dealing with strongly subsampled data, our first aim is to get a solution that fits the projection data as much as possible. This suggests turning to CS reconstruction that targets an unbiased fit to the subsampled data. We thus carried out CS reconstructions relying on a homotopy strategy with soft background subtraction, TV minimization and thresholding in the wavelet domain. Reconstruction settings are indicated in Tab. Experiment Algorithm CS hyperparameter while preserving high frequency information from FBP reconstruction. However, we notice that in the region of air, at an intensity level inferior to the one of the soft-tissues, we get a reconstruction that is similar to standard ℓ 2 reconstruction. We also notice that the cerebellum appears not as resolved that in FBP reconstruction. CS-tv corrected for streak artifacts and recovered the bony structures better than CS-sbs. However, despite the use of a homotopy strategy, it still resulted in a patchy appearance which would not be acceptable for the clinical practice.

n • 1 iFBP+ n • 2 CS-sbs λ (s) sbs = 3000 to 0 n • 3 CS-tv λ (s) tv = 1000 to 0 n • 4 CS-w λ (s)
We monitored the convergence for both CS-sbs and CS-tv. Curves of log d r with respect to the number of iterations over the whole image are plotted in Fig. 4.12. CS-sbs and CS-tv achieved a discrepancy level that is below the one of iFBP+, whose convergence is shown as a black dashed line. It reflects the successful mitigation of streak artifacts by both sparsity constraints.

Extension to non-ℓ 1 -norm constraints

Let us now replace TV minimization with anisotropic diffusion in a proximal splitting scheme with iFBP and homotopy. Figure 4.13 shows the reconstructed image. The reconstruction is very close to CS reconstruction using TV: high contrasts are preserved, while low contrasts tend to form piece-wise constant regions. Small differences can be observed, although it is difficult to decide which operator between TV and anisotropic diffusion yields the best reconstruction result. This result is particularly interesting since it demonstrates that a reconstruction strategy that neither uses random measurements nor relies on the minimization of a ℓ 1 -norm can provide to some extent the expected benefits of CS [START_REF] Langet | Nonlinear diffusion constraints for reconstructing subsampled rotational angiography data[END_REF].

Uniform subsampling with clinical data

In this study, we investigate the tomographic reconstruction from 150 views only, since it is the least favourable sampling case that can be met in clinical practice. Clinical data are acquired on an Innova 4100 C-arm system (GE Healthcare, Chalfont St. Giles, UK). The system records equiangular cone-beam 2D views at 30 frames/s during an approximately 200 • rotation at 40 • /s delivering 150 views in total. Reconstruction settings are indicated in Tab. 4.5. 4.5: Generic parameters of iFDK-derived algorithms.

Subtracted data

The first dataset we study is a DSRA exam of cerebral vessels. A subtracted scan was provided by subtraction of the mask scan from the contrast scan. Figure 4.14 displays a slice from reconstruction of the full subtracted scan, in which sparse contrast-enhanced vessels and non-sparse perfused tissues (grey areas surrounding the vessels) also referred to as parenchyma are visible. In this example, one sees that the subtracted volume does not contain sparse areas only, which shows the importance of an unbiased reconstruction. FDK reconstruction of Fig. 4.14a presented streak artifacts that are issued from the vessels and are more prominent at the periphery. Both the positivity constraint (Fig. 4.14b) and SBS (Fig. 4.14c) promoted sparsity of the reconstructed structures, thus reducing streak artifacts with respect to FDK reconstruction. Roughly speaking, two steps can be identified in the streak removal mechanism that is implemented in iFBP+:

• the positivity constraint removes negative streaks;

• while iFDK back-forward projection step removes the complementary positive streaks. On the other hand, visual analysis of the respective slices showed that iFDK+ reconstruction enhanced high-frequency information, which yields a noisy visual appearance. This effect disappeared using CS -sbs reconstruction, that then compares favourably with FDK reconstruction. To confirm this visual appreciation, we computed the mean and standard deviation in a region of interest of 930 voxels within these perfused tissues and found 116 ±142 HU in Fig. 4.14b for positivity and 96 ±55 HU in Fig. 4.14c for SBS, a twofold Signal-to-Noise-Ratio (SNR) increase for CS -sbs over iFDK+ (1.8 vs. 0.8). Comparison of profile curves drawn in Fig. 4.15 through vessels revealed no differences in terms of resolution. It is however obvious that less intensity variations in the background structures were introduced by using CS -sbs. Hence, soft background subtraction (SBS) provided to a certain extent the expected benefits of the CS theory: reduced streaks. Such gains would otherwise be achievable only by increasing the number of projections, which was not possible in these clinical settings, rather than increasing the X-Ray dose of each projection image.

Non subtracted data

In a second study we consider a clinical dataset which is a scan of about 150 views of slightly opacified carotid arteries. The presence of a tooth metallic implant in the imaged data results in strong horizontal streaks. Soft tissues show small streaks coming from the bones.

Cone-beam subsampling

In the previous clinical case, streaks were corrected despite the low number of views. This implicitly assume that the streaks would not have been visible with an acquisition at a higher frame rate. In the following, we consider a clinical non-injected acquisition of 600 views of the head and provide CS reconstruction with homotopy of those 600 views as well as subsets of the views so as to recreate uniform subsampling. Figure 4.19 demonstrates that two types of streaks remain at 600 views:

• so-called metal artifacts which correspond to data that are corrupted by very high noise due to the density metal,

• CB artifacts since the circular orbit does not fully sample the objects whatever the number of views.

Only the CB artifacts get corrected by CS. In addition, it is shown in Fig. 4.20 that streaks originated from an object outside the FoV, that appear for low sampling rates, cannot get corrected by CS. In Fig. 4.21 and Fig. 4.22, the CB artifact (horizontal black streak) is always corrected by CS whatever the number of views. The correction obtained with the TV penalty is however superior to the one with SBS at very low number of views. In particular, we observe that sinuses are blurred with SBS (at the air-tissue interface), while their resolution is preserved with TV as we observed on the simulated data. This can also be seen in Fig. 4.23 where reconstruction with TV yields, in particular, better depiction of the alveolar structures of the temporal bone than SBS.

Limited range subsampling with clinical data

We now consider the case of data that are corrupted by contrast variations in the injected vessels, the acquisition pattern associated to a phase is reduced to a limited angle. As discussed in Sec. 2.1.3, FDK reconstruction of limited range data is strongly artifacted and should be avoided by reconstructing a full range of inconsistent data instead. Here we challenge this strategy by applying the ℓ 1 -penalized iFDK with homotopy to projection data with limited angular coverage.

Subtracted data

We consider the same DSRA exam of cerebral vessels that was used for producing the reconstruction results presented in Fig. 4.14. During the acquisition of the contrast scan, the right vertebral artery (RVA) did not appear opacified during the first half of the rotation.

Non subtracted data

Application on non-subtracted data is shown in Fig. 4.26. The same slice as in Fig. 4.24 was reconstructed without subtraction. Besides what we already noted in the subtracted case (Fig. 4.24(a)), the nonsubtracted FDK reconstruction (Fig. 4.26(a)) showed that the horizontal and vertical black streaks associated to the RVA hid the underlying tissues. Limitating FDK reconstruction to the second half of the scan (Fig. 4.26(b)) yielded no visual improvement of the RVA while degrading both LVA and background. On the contrary CS-sbs (Fig. 4.26(c)) recovered both arteries, as in the subtracted case. None of these solutions is entirely satisfactory because they do not recover the non-sparse background which can only be reconstructed from the entire scan.

Discussion

We proposed in this study algorithms that enable the generation of reconstruction results mitigating cone beam aliasing artifacts with a relatively low computational cost. Another key aspect in tomographic reconstruction is noise management. To take noise into account, the computation of the fidelity term has often been handled differently to the least-square approach:

• multiplicative ART (MART) relies on an entropic criterion [START_REF] Gordon | Algebraic Reconstruction Techniques (ART) for threedimensional electron microscopy and x-ray photography[END_REF];

• statistical reconstructions with Poisson laws are based on the maximization of the log-posterior probability (maximum a posteriori or MAP) among which the expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Shepp | Maximum Likelihood Reconstruction for Emission Tomography[END_REF][START_REF] De Man | An iterative maximum-likelihood polychromatic algorithm for CT[END_REF].

The validity of CS is established with the least square data fit only. The generation of sparse solutions with other data fit criterion has been explored in [START_REF] Ting | Sparse Image Reconstruction for Molecular Imaging[END_REF][START_REF] Dobigeon | Myopic sparse image reconstruction with application to MRFM[END_REF][START_REF] Pustelnik | Méthodes proximales pour la résolution de problèmes inverses: application à la tomographie par émission de positrons[END_REF] and is beyond the scope of this work.

Our approach is applicable to any tomographic dataset to provide to a certain extent the expected benefits of compressed sensing theory: reduced streaks. Such gains would otherwise be achievable only by increasing the number of projections, which is not possible in these clinical settings, rather than increasing the X-Ray dose of each projection image. The linear decrease of the soft background subtraction threshold provided superior reconstruction of the non-sparse background of sparse high-intensity structures, resulting in more uniform tissue depiction when affected by streaks from either teeth, bones or opacified vessels, at equal resolution. It has been shown to efficiently mitigate subsampling streaks due to high-intense structures, but it is less efficient when dealing with artifacts originating from soft tissues. In addition, it is not symmetrical: it preserves positive peaks (i.e. the contrast-enhanced vessels) but not the negative ones (air-soft tissue interface). The best reconstruction was obtained using TV with homotopy. Total variation is symmetrical and has been proven to remove streak artifacts efficiently, but at the expense of a patchy look. It is better than SBS for very low number of views. Note that streaks originating from metal or objects outside the FoV cannot be removed.

We used a very simple approach for the thresholding of wavelet coefficients and observed that reconstruction results are corrupted by compression artifacts. When the sampling is more pronounced, thresholing of the wavelet transform is not able to deal with severe streaking artifacts: in particular, it captures the vessels, but is not able to discriminate streaks from other structures (while TV does), and hence is not applicable to less favourable subsampled problem. However, we are aware this implementation is surely not optimal and more complex approach such as translation invariant wavelet frame [START_REF] Mallat | A review of bandlet methods for geometrical image representation[END_REF] or tree frame1 [71] could be considered for some improvements.

With limited range, none of the proposed CS approaches recovers the background, only a full 180 • reconstruction can restore it. This is the subject of the next chapter.

Chapter 5

Reconstruction of dynamic data

The generic reconstruction framework that we discussed in the previous section finds an important application in dynamic imaging. We formulated in Sec. 2.1.3.3 the problem of dynamic data acquisition with C-arm systems as a time-sequential sampling problem, where only one view can be taken at a time and the time interval between successive views is related to their angular interval. The reconstruction artifacts that arise due to this insufficient acquisition rate are then clearly identified as subsampling artifacts. This chapter aims at extending the proposed reconstruction framework for the correction of motion artifacts in rotational angiography, assuming that the support of the temporal variations is sparse in both the spatial and the temporal domains. Redundancy within time series is exploited either for imaging moving objects or for imaging objects whose contrast is rapidly changing, such as in dynamic contrast uptake studies, and contrast angiography.

Dynamic problem statement

To begin this chapter, let us shortly introduce the notations that are useful for describing the reconstruction of dynamic data. Let us describe a dynamic object by set of volumes f dyn ∈ R K×M , where volume f m ∈ R K refers to a given phase of the object among M possible phases: within f dyn during the acquisition correspond to a set of phases f m . The key point in this context of C-arm imaging is that the time series is sampled by a unique acquisition mode of N angular positions, for instance in the 2-phase case illustrated in Fig. 5.1:

f dyn =    f 0 . . . f M -1    . Coefficient (f m ) k is
t = {f 0 , f 0 , f 0 , f 0 , f 1 , f 1 , f 0 , f 0 , f 1 , f 1 , • • • } .
We split measurement vector p into M subsets so that vector p m ∈ R J×Nm contains the N m projections at the angular positions assigned to phase f m only. For simplicity of notation, we define measurement vector p dyn such that:

p dyn =    p 0 . . . p M -1    .
Let us then denote R m ∈ R J×Nm × R K the subset of projection matrices from R relating phase f m to subset p m (i.e. such that p m = R m f m ) and R dyn the block-diagonal matrix such that:

R dyn =    R 0 . . . R M -1    .
With these notations, the dynamic reconstruction problem can be written:

p m = R m f m p dyn = R dyn f dyn . (5.1)
Phase-wise FBP reconstruction fdyn,FBP ∈ R K×M is obtained by phase backprojection R T dyn of the filtered projections: fdyn,FBP = R T dyn Dp dyn , (

or minimization of the least-square criterion:

arg min

f dyn ∈R K×M 1 2 R dyn f dyn -p dyn 2 2,D = Q dyn (f dyn ) ≡ M -1 m=0 arg min fm∈R K 1 2 R m f m -p m 2 2,D = Q m (f m ) (5.3)
while static FBP reconstruction fstatic,FBP ∈ R K corresponds to the weighted sum of the phase-wise FBP reconstructions:

fstatic,FBP = R T Dp = M -1 m=0 w m fm,FBP , (5.4)
where

w m = π N • Card {E m } as defined in Sec. 2.1.3.

Motion in rotational angiography

In this section, we give an overview of the reconstruction strategies that have been previously proposed in the literature for the management of motion in rotational angiography.

Motion compensation

It is possible to reduce the number of unknowns of the linear system given in Eq. (5.1) by modelling the motion between phases. Let us denote T m : R K → R K the time-dependent transformation that models the deformation of phase f r that is arbitrarily chosen as reference into phase f m and T -1 m its inverse. With these notations, let us denote T : R K → R K×M the transformation operator such that:

f dyn = T f r =    T 0 f r . . . T M -1 f r    .
Parametric models are employed to address global deformations with a small number of degrees of freedom. To compensate for rigid body motion, Lin et al. proposed in [START_REF] Lin | Compensation for movement in computed tomography equipment[END_REF] an iterative method that estimates the centroid of the imaged object from CB projections. To compensate for the respiratory motion, Lu et al. [START_REF] Lu | Tomographic motion detection and correction directly in sinogram space[END_REF] proposed to rescale the projection data by means of scaling parameters that are derived by tracking markers or fiducial points and then apply standard FBP, while Crawford et al. [START_REF] Crawford | Respiratory compensation in projection imaging using a magnification and displacement model[END_REF] rather developed a dedicated algorithm that incorporates a magnification model into the reconstruction.

To address local displacements of the vessels that are deformable structures, non-linear (dense) models are favoured. In particular, Blondel et al. [START_REF] Blondel | Reconstruction of coronary arteries from a single rotational X-ray projection sequence[END_REF] proposed a motion-compensated reconstruction of the coronary arteries for C-arm CT. The approach consists in generating a 3D model of the artery centrelines, from which a 4D motion vector field is estimated; the vector field is then used to modify the projection matrices. It has also been successfully applied to the reconstruction of coronary stents [START_REF] Perrenot | Motion correction for coronary stent reconstruction from rotational x-ray projection sequences[END_REF]. The estimation of the 4D models relies on the detection of the arteries within the 2D projections, that are then tracked over the acquisition. This method cannot be generalized to other anatomies where the vasculature is either too complex, as in the brain, or too attenuated by the surrounding structures as in the liver. We have seen in Sec. 2.3.4.3 that it is generally not possible to incorporate motion models into the analytical reconstruction formulas. Reconstruction strategies rather handle motion compensation than phase-wise reconstruction post-processing. Let us assume that local correction is a valid approximation in tomography, that is to say motion compensation can be incorporated within the backprojection step by shifting the voxel to be reconstructed according to the motion displacement field. It follows that motion-compensated reconstruction is the weighted sum of the deformed phase-wise FBP reconstructions [START_REF] Schäfer | Motion-compensated and gated cone beam filtered back-projection for 3-D rotational X-ray angiography[END_REF][START_REF] Li | Motion correction for improved target localization with on-board cone-beam computed tomography[END_REF]:

fstatic,FBP+T = fr,FBP+T ≃ M -1 m=0 w m T -1 m fm,FBP . (5.5)
Such a local correction does not handle the global artifacts that remains in other regions [START_REF] Li | Motion correction for improved target localization with on-board cone-beam computed tomography[END_REF], i.e. regions that are static but still suffer from the artifacts that are generated by the moving structures. An iterative strategy is thus preferred to incorporate the motion model, that yields a reconstruction with superior motion artefacts correction [START_REF] Li | Improved iterative algorithm for sparse object reconstruction and its performance evaluation with micro-CT data[END_REF].

The level of complexity of current motion models is however not adequate for addressing standard C-arm imaging marred by accidental motion artifacts or contrast flow variations. More importantly, as illustrated in Fig. 5.2, the perfect knowledge of the motion does not necessarily induce sufficient angular sampling of the object and hence, does not guarantee the feasibility of the reconstruction. In this toy example, the 2D phantom of Fig. 5.2(a) is simulated in rotation at the same speed as the system gantry acquiring measurements in parallel geometry. The FBP reconstruction from the motion corrupted projection data set of Fig. 5.2(b) yields the artifacted volume shown in Fig. 5.2(c), while the motion-compensated reconstruction is the FBP reconstruction of a single sinogram profile as illustrated in Fig. 5.2(d). It is obvious that perfect knowledge of the motion does not improve the angular range and hence, does not enable the tomographic reconstruction of the original phantom.

Motion reconstruction

Compressed sensing reconstructions have been investigated for dynamic reconstruction in cardiac imaging [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF][START_REF]Projection-based motion compensation for gated coronary artery reconstruction from rotational x-ray angiograms[END_REF] and radiotherapy [START_REF] Bergner | A comparison of 4D cone-beam CT algorithms for slowly rotating scanners[END_REF]. We have shown in Sec. 4.7 that, except for the particular case of subtracted data, the CS reconstruction of uniformly subsampled volumes does not recover the background structures. Consequently, a CS phase-wise strategy solely relying on spatial penalties generally does not allow for reconstructing dynamic data. Sparsity must be expressed in terms of both spatial and temporal penalties. The temporal sparsity is illustrated on the angiographic sequence that we presented in Sec. 3.4.2.2, of which we display a zoomed detail of the position of the catheter for two different time-points in Fig. 5.8(a), as well as the sinogram of the difference in Fig. 5.8(b). Temporal sparsity, that is applicable even in absence of a motion model, was explored with CS, yielding a sparse approximation to the solution through spatial constraints (typically wavelets or total variation filtering) and temporal redundancy [START_REF] Sidky | A constrained, total-variation minimization algorithm for low-intensity x-ray CT[END_REF][START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF][START_REF] Jia | 4D Computed Tomography Reconstruction from Few-Projection Data via Temporal Non-local Regularization[END_REF] and temporal sparsity in the Fourier domain considering periodic variations for cardiac MRI [START_REF]Sparse MRI[END_REF].

Because the temporal variations within the anatomy are localized, the motion and contrast variations of the contrast-enhanced vessels in rotational angiography can be described by means of a two-layer model:

• a static background that is common to all phases f m ;

• a dynamic layer that contains structures in motion and occluding the background. The support of the moving structures is specific to each phase f m . This layer is sparse.

The collected views thus correspond to the projection of the moving structures superimposed on the projection of the static background. The key idea consists in identifying and differently reconstructing the voxels of each layer:

• voxels that belong to the static background are reconstructed from all the measurements (i.e. the full scan);

• while voxels that belong to the dynamic layer are reconstructed from measurements of a single phase.

Two cases can be considered:

• the layers can be separated in the projection domain: given a segmentation of the moving objects in the 2D views, it would be possible to suppress the projection associated to these objects and to perform a reconstruction of the background only. This strategy requires explicit vessel segmentation;

• the layers superimposition is lost in the projection domain and the redundancy of the static structures is used instead: provided an estimation of the background, it would be possible to detect the moving structures by subtracting the reprojection of the background from the projection data.

Redundancy of the background structures is similar to video compression (e.g. MPEG-1), that uses key frames of a video sequence to encode for the (sparse) differences of all frames with respect those key frames, possibly relying on motion estimation to reduce sparsity even further.

Layer separation in the projection domain

Non-linear identification of the support of the moving structures directly in the 2D views was proposed in the auto-adaptative phase-correlated (AAPC) reconstruction that was developed by Bergner et al. [START_REF]Autoadaptative phase-correlated (AAPC) reconstruction for 4D CBCT[END_REF] in the clinical context of radiotherapy. In this approach, the motion field between two projections is computed through an optical flow algorithm. Optical flow techniques rely on a local intensity conservation postulate (i.e. assumption that intensity remains constant from one view to the other), an assumption that is not valid with (projective) X-ray images. The optical flow algorithm used in AAPC relies on a global intensity conservation postulate instead (i.e. assumption that the sum of all intensities remains constant from one view to the other). However its use is limited to non-truncated data. In the case of truncated data, the measurements also include attenuations from objects out of the imaged FoV, that may not be the same for all angles. The velocity change of direction is detected so that a pixel is said static if its velocity has a direction quasi-constant in a given observation window (to take the rotation of the acquisition system into account). In this way, each pixel having a velocity with a non-null component along the acquisition system rotation axis is considered in motion. For more robustness with respect to noise, only the most important magnitudes are kept. It is a suitable criterion for cyclically moving objects. Based on this segmentation, each pixel in a view is estimated as being either motionless or motion-affected and is assigned a weight. Pixels that are assigned too low a weight are dynamically compensated by interpolating the gaps within the backprojection integral with trapezoids using the high weight values. Standard FBP reconstruction is finally performed.

Redundancy of the background

Generally we do not have an external estimate of the background. Therefore we have to estimate the background directly from the data set. In particular, in the case of sequences with numerous and strong artifacts, obtaining a reference image is as difficult as segmenting the sequence, in particular in the case of noisy sequences. A simple idea is to compute a prior image, which has been investigated in the following works.

Mc Kinnon-Bates (MKB) An iterative method that aims at mitigating motion artifacts is the twostep algorithm algorithm that was proposed by McKinnon and Bates [START_REF] Mc Kinnon | Towards imaging the beating heart usefully with a conventional CT scanner[END_REF], referred to as MKB. The static FBP reconstruction fFBP is produced, then forward-projected at the same angles as those of the original projection data in order to generate a set of views that is contaminated with the motion artifacts.

The difference due to the motion is enhanced by subtraction of the computed views from the original projection data, then reconstructed for each phase using phase-wise FBP. The difference images are added to static FBP to generate the final phase-wise reconstruction:

fm,MKB = fstatic,FBP + R T m D(p m -R m fm,FBP ) (5.6)
In the static regions, the final phase-correlated images differ little from static FBP (nominal zero values), while in regions where motion is present, the final reconstruction more closely resembles the phase-wise FBP. One limitation of the MKB algorithm is that, although it provides an efficient means of mitigating the streaks that are caused by static structures, it cannot remove the streak artifacts of moving structures so that corrected cardiac CT is not clinically acceptable. Another limitation of the MKB algorithm is that the noise in the corrected images is then primarily determined by the subsampled phases. Thus, the final image has a lower contrast-to-noise ratio (CNR). Some refinements for noise suppression have been proposed by Garden et al. [START_REF] Garden | 3-D Reconstruction of the Heart from Few Projctions: A Practical Implementation of the McKinnon-Bates Algorithm[END_REF]. A final limitation is the presence of truncation artifacts due to the forward-projection operation for which Zheng et al. proposed data interpolation [START_REF] Zheng | Fast 4D Cone-Beam CT Reconstruction Using the McKinnon-Bates Algorithm with Truncation Correction and Nonlinear Filtering[END_REF]. The McKinnon-Bates (MKB) algorithm was successfully used to correct for streaking artifacts due to the respiratory motion in lung data for image-guided radiation therapy [START_REF] Leng | Streaking Artifacts Reduction in Four-dimensional Cone-beam Computed Tomography[END_REF].

Highly constrained backprojection (HYPR) Mistretta et al. [START_REF] Mistretta | Highly constrained backprojection for time-resolved MRI[END_REF] proposed for DSA or digital subtraction rotational angiography (DSRA) an empirical approach, called HighlY constrained back-PRojection (HYPR). It uses a so-called 'composite image' that is actually a static FBP reconstruction fstatic,FBP to constrain the frame-wise reconstructions. To this aim, the composite image is forwardprojected, then backprojected. Backprojection of the projection is related to a given phase. Note that there is no use of the ramp filter so as not to enhance high frequencies that would make the weighting explodes. Each view is assigned a weighted image that is defined as the ratio projection over projection of the composite image. The motion compensated image is defined as the product of the prior image by the weighting map. HYPR is claimed to work for the reconstruction of vascular structures, but lacks rigorous mathematical background.

Prior image constrained compressed sensing (PICCS) Formalizing the FBP static reconstruction used in MKB to the more general idea of a prior image -estimation of the static background structures-, jointly with ℓ 1 -norm minimization that allows for discriminating the static support from the varying support, the Prior Image Constrained Compressed Sensing (PICCS) algorithm has been proposed by Chen et al. [START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF][START_REF] Theriault Lauzier | Prior image constrained compressed sensing: Implementation and performance evaluation[END_REF][START_REF] Chen | Time-resolved interventional cardiac C-arm cone-beam CT: An application of the PICCS algorithm[END_REF] to pool together inherently incoherent information, and recover the dynamic variations that affected the data. This approach is particularly seducing since it relies on using a temporal constraint based on a prior image (that may be the static reconstruction), so that the optimization process, schematically, determines motionless areas from the complete data set, while motion-blurred areas are determined from subsampled reconstruction similar to what is done in MKB. With respect to MKB, PICCS allows for the easy incorporation of additional prior information into the optimization process, and hence the proposed approach includes total variation minimization for mitigating appearance of subsampling streak artifacts in the reconstruction results. A comparative study demonstrated significantly higher quality using the PICCS algorithm [START_REF] Bergner | A comparison of 4D cone-beam CT algorithms for slowly rotating scanners[END_REF]. Though this idea of comparing a reference image with the current image is very intuitive, it is not always clinically applicable. In particular, we have seen in Sec. 3.2 that with intraarterial contrast injection, motion blur is not the only degradation present in the standard FBP static reconstruction: any part of the image may be degraded by intense streaks due to the presence of inconsistent high-intensity vessel projections in the data, making the static reconstruction a poor estimation of the background structures. In this work, we shall propose a generalization of the PICCS algorithm that is adapted to rotational angiography with intraarterial injection.

Compressed sensing for digital subtraction rotational angiography

We have seen that CS enables a more efficient use of the collected samples. It also suggests to redesign the acquisition protocol. Let us explore this aspect of CS in the particular case of DSRA.

Acquisition protocols

Let us consider a protocol that collects N equiangularly spaced views for the contrast and mask volumes.

Homogeneous protocol To reconstruct a satisfying subtracted volume, the mask and contrast scans are currently collected with the homogeneous protocol that is schematized in Fig. 5.3(b). In this protocol, the trajectories of the mask and the contrast scans sample the same set of angular positions:

Φ C = Φ M = {φ n = n • ∆φ | n = 0 • • • N -1} (5.7)
The acquisition with identical parameters allows for the straightforward removal of redundant background structures and their associated streaks as shown in Fig. 5.4(a). Filtered backprojections of the mask and the contrast suffer from the same undersampling artifacts and the same noise level. They cannot be combined a posteriori to either reduce noise or artifacts because vessel streaks would then propagate into the mask volume.

Heterogeneous protocol Langet et al. [START_REF] Langet | Sparsity Constraints and Dedicated Acquisition Protocols for Improved Digital Subtraction Rotational Angiography[END_REF] proposed a novel acquisition protocol for DSRA that is schematized in Fig. 5.3(b) and in which each scan samples two sets of interleaved angular positions defined by:

Φ M = {φ n = n • ∆φ | n = 0 • • • N -1} Φ C = {φ n = n • ∆φ + δ | n = 0 • • • N -1} (5.8)
where δ =]0, ∆φ[ is the angular shift. Acquiring heterogeneous scans allows for increasing angular sampling and therefore image quality of the non-opacified structures through a simple average operation (since they are sampled more often). In particular, with δ = 0.5∆φ, background structures are sampled with 2N equiangular-spaced views, whereas the vessels are sampled with N equiangular-spaced views as illustrated in Fig. 5.5(c) with respect to Fig. 5.5(a) and Fig. 5.5(b). However, as can be seen in in Fig. 5.5(c), vessel streaks cannot be removed. In addition, when subtracting the two volumes, the streak artifacts corresponding to the nonopacified structures are amplified as shown in Fig. 5.4(b), therefore considerably altering the reconstruction of the subtracted volume.

Alternative protocols Sampling improvement could be achieved by acquiring the mask with twice as many frames as the opacified. Mask projections are thus acquired at all opacified angles and at all intermediate angles:

   Φ M = φ n = n • ∆φ 2 | n = 0 • • • 2N -1 Φ C = {φ n = n • ∆φ | n = 0 • • • N -1}
This however would increase the total dose by 50%. On the other hand, such a method seems to offer the possibility of 2D subtraction since there would exist a mask image for each opacified image. However, C-arm CT systems are operated in 3D at their maximum frame rate. Therefore, in order to double the number of frames of the mask scan, the mask rotation speed must be divided by two. If mask and opacified rotation speeds differ, there is no guaranty that it will be possible to acquire a mask image at the exact same positions as the opacified one. Therefore, 2D subtraction is not guaranteed in practice.

The same goals could be achieved through the use of a prior examination obtained from a separate protocol possibly with another modality such as CT or MRI, from which a high-resolution mask can be computed. The fundamental difference with our approach is that registration is required when using a prior examination, which is a difficult problem in itself, whereas we assume here acquisition with a single protocol fast enough so that the patient can be maintained still so that registration is not necessary.

f M f C Scan start Scan end f M f C Scan start Scan end (a) (b)

Two-phase CS reconstruction algorithm

Clearly DSRA is a particular case of dynamic rotational angiography, in which the collected data encode for two phases:

• the mask phase f M ∈ R K that is imaged without vessel enhancement;

• the contrast phase f C ∈ R K that is imaged while the vascular structures in the FoV are fully opacified.

Let us denote f S ∈ R K the subtracted volume that is obtained by subtraction of the mask from the contrast volume:

f S = f C -f M .
Using this dynamic interpretation, we define:

f dsra = f C f M ,
and the reconstruction problem now fits in Eq. ( 5.1) with standard iFBP reconstruction computed in Eq. (5.3). Temporal penalty Assuming that the contrast-enhanced vessels that are contained within the subtracted volume f S are naturally sparse structures (see Sec. 4.4), we address the DSRA reconstruction by promoting sparsity of f S or, symmetrically, redundancy of the nonopacified structures captured in both the mask and the contrast scans, whatever the acquisition protocol. To mix background information from the contrast and the mask without losing vessel quantification, we define an orthonormal temporal transform H t such that f S is an isolated component of product H t f . Hence, H t is defined as the operator associated with the 1D Haar wavelet transform:

H t = 1 2 1 2 1 -1
It follows:

H t f dsra = f C + f M 2 f S .
Subtracted volume f S is captured by the detail subband of the Haar wavelet transform. Thus, CS t reconstruction can be handled as a joint iterative reconstruction of the contrast and the mask volume, while maintaining a sparse constraint of their coupling through H t . Moreover, positivity of all three volumes f M , f C and f S is ensured.

fdsra,χt λt = arg min

f dsra ∈R K×2 Q dsra (f dsra ) + ι + (f dsra ) + λ t f S 1 + ι + (f S ) = λ t • χ t (f dsra ) , (5.9)
where Q dsra (f dsra ) is the data fidelity term as in Eq. ( 5.3). As presented in Sec. 4.2.2, the proximal operator of χ t (f S ) is soft background subtraction. Since the sparsity assumption on the subtracted volume is valid, a regularization strategy is proposed and λ t is assigned a fixed value. Reconstructing a sparse subtracted volume implies it cannot contain streaks. As a consequence, these streaks are also removed from the mask and opacified volumes. Such an approach eliminates the need for identical sampling of the mask and the contrast scans.

Spatial penalty

The maximal sampling that can be achieved for f C and f M cannot exceed the combination of the mask and the contrast sampling, which may still result in a subsampled problem. To remove streaks from the contrast volume and prevent the transfer of vessel-related information (including streaks) into the mask, since both volumes are not sparse, we penalize both the mask and the contrast volumes ℓ 1 -norms using a spatial sparsity constraint. Here we propose to combine the temporal sparsity constraint with soft background subtraction, and thus CS sbs+t reconstruction problem can be constrained with the following penalty: fdsra,χ sbs +χt λ sbs ,λt = arg min

f dsra ∈R K×2 {Q dsra (f dsra ) + λ sbs • χ sbs (f dsra ) + +λ t • χ t (f dsra )} (5.10) A set of decreasing λ (s)
sbs values, as described in Sec. 4.36, is used. Note that any proximal operator prox χ A that promotes spatial sparsity of the data can replace SBS (e.g. TV minimization with A = ∇):

fdsra,χA+χt λ A ,λt = arg min f dsra ∈R K×2 {Q dsra (f dsra ) + λ A • χ A (f dsra ) + λ t • χ t (f dsra )}
(5.11)

Reconstruction results

The following investigates image quality of subtracted rotational angiography by using heterogeneous scans with adequate sparsity-based reconstruction.

Study on Simulated Data For our simulation, we used a 512 × 512 cerebral CT cross-section as mask image. We simulated the contrast image by adding to the mask synthetic disks that represent opacification. The values vary from 2000 to 3000 HU for the simulated injected vessels and around 1200 HU for the simulated parenchyma, while soft tissue values (around 1000 HU) and bone values (around 2000 HU) are those of the original CT slice. We produced interleaved mask and contrast scans in parallel geometry from 150 integral lines that are equiangularly spaced over 180 • with angular shift δ = 0.6 • for the contrast scan. We compare the reconstruction quality of the background structures in Fig. 5.6. Standard reconstruction of the mask with a double sampling (300 views) is shown in Fig. 5.6a as reference: its streak level is the lowest level that can be achieved with the approach we developed (d r = 9 HU). Subsampled standard reconstruction (150 views) yields streak artifacts that, in particular, makes cerebral sulci visualization difficult, as shown in Fig. 5.6b. This degradation was confirmed by RMSD value d r = 25 HU (see Sec. 3.1.4.3 for RMSD definition). The contrast image (see Fig. 5.6c) presents additional streaks due to the injected vessels. Figure 5.6d displays mask reconstruction penalized by the temporal constraint χ t (f dsra ) only with weight λ t = 5. It allows for removing the background streaks and for recovering a background resolution similar to the reference one. Nevertheless, we notice that a small amount of the highest contrasted structures (i.e. the injected vessels) is transferred from the contrast to the mask, and so do the very thin associated vessel streaks, which affect the whole image (d r = 15 HU) and are well visible near the skull bone on the right of the image. Such a transfer is not observed for lower opacification contrasts (e.g. parenchyma). Vessel marks and streaks are not visible anymore in Fig. 5.6e that was produced with CS sbs+t while using a spatiotemporal penalty. Streak level was d r = 11 HU, which confirms the effective image quality improvement. True vessel intensity is recovered, which would not be the case if we simply computed the average reconstruction of the mask and the contrast volume. For all reconstructions, sparsity of the subtracted volume (not shown here) is preserved by construction: it perfectly fits the sparsity model, an assumption which will no longer hold for real data.

Study on Clinical Data

In DSRA routine, the C-arm system records projections at 30 frames/s during an approximately 200 • rotation at 40 • /s delivering 150 views in total. Since there exists no clinical protocols with interleaved scans, we built an interleaved acquisition pattern with 75 views for each scan by taking one view every two views with a shift of one view when starting the contrast scan as illustrated in Fig. 5.3(b). Figure 5.7 compares standard reconstruction (iFDK+, displayed in the first row) and CS sbs+t (displayed in the second row). Similarly to previous simulations, our algorithm yields significantly fewer streaks than standard reconstruction at half sampling in clinical data as well. The resulting resolution improvement is best seen in the petrous part of the left temporal bone (right side of the slice, detail zoomed in Fig. 5.7b) that contains the inner ear: thin details such as tympanic cavity, canals, and sutures are more accurate. Computation of the mean and standard deviation in a region of interest of 900 voxels within soft tissues confirms these observations: 1308 ±243 HU is found for standard reconstruction, while 1077 ±166 HU is found for our reconstruction, which corresponds to a 20% Signal-to-Noise-Ratio (SNR) increase (6.5 vs. 5.4). Note that 1116 ±112 HU corresponding to a SNR value of 10.0 is found for standard reconstruction with 150 views (not shown here). Looking at the subtracted volume shown in Fig. 5.7c, we notice that even if our approach yields sparser structures than standard reconstruction at half sampling, the whole volume is not sparse, which limits the background mix between the mask and the contrast volumes.

Discussion

The results that were obtained on both simulated and clinical data showed that our approach outperforms standard reconstruction in terms of background restoration and streak removal. Moreover, this quality improvement is performed at a reasonable cost of 20 iterations and a small change in the acquisition protocol so that both scans sample interleaved angular positions. Vessels reconstruction with a prior of sparsity is not a new idea, and it is being used for reconstructing sparse structures. Our strategy takes advantage of a sparse prior, despite the fact that the two scans corresponds to different objects and that neither the mask nor the opacified volumes are sparse, but only their difference, to improve the 3D image quality of the mask, opacified and subtracted volumes. The improvements in image quality have been investigated for subtracted cerebral data that contain vessels only. It would be more challenging to apply the CS approach to abdominal data, since motion is more likely to occur between both acquisitions. The validation of the CS approach in the case of subtraction that contains parenchyma as well is still a work-in-progress, that requires more validation. This example is simple, since two phases are considered only, there is no overlap between phases, no error in phase identification and there are two rotations instead of one only. As we shall see in the following, motion reconstruction is more challenging.

Compressed sensing for motion reconstruction

We have seen that, with intraarterial contrast injection, any part of the static FBP reconstruction may be degraded by intense streaks due to the presence of inconsistent high-intensity vessel projections in the data, making the static reconstruction a poor prior image in PICCS algorithm proposed by Chen et al. [START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF]. CS-sbs reconstruction that we detailed in Sec. 4.4.3 could mitigate subsampling artifacts while preserving the overall aspect of FDK reconstructions, a point of importance for the clinical practice. However, this strategy is not applicable to removing streaks due to incoherent static data and thus cannot provide an improved prior image. We therefore adapt our framework to the reconstruction of a series of volumes based on combining the SBS spatial constraint to a temporal constraint that enforces the sparsity of the difference between time points. This contribution is applied to two types of temporal inconsistencies within injected vessels of the abdominal regions:

• small displacements induced by the blood flow pulsatility,

• and contrast flow variability during the scan due to delayed opacification.

Multiphase CS reconstruction algorithm

CS-sbs enables the independent reconstruction of each phase so that the sampling of the static background structures is significantly reduced. To recover full sampling of the background we have to take inter-phase correlations into account: fdyn,χA+χt

λ A ,λt = arg min f dyn ∈R K×M Q dyn (f dyn ) + λ A • χ A (f dyn ) + λ t • χ t (f dyn ) , (5.12)
where the penalization term λ t •χ t (f dyn ) describes the temporal sparsity of the data. Spatial sparsity can be enforced through soft-background subtraction or total variation minimization as previously presented.

In PICCS [START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF], temporal sparsity is handled by assuming sparsity of the difference to a prior image f P . Such an approach can be nested within a proximal splitting algorithm with a penalty that is defined as: We denote CS-piccs the reconstruction algorithm that uses a prior image chosen equal to the static FBP reconstruction as in [START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF]: f P = fstatic,FBP , and SBS. Note that with this definition the prior image corresponds to the weighted mean calculated from all time-point FBP reconstructions as given in Eq. (5.4). Using proximal operators, we can use a larger class of filtering penalties to solve the dynamic problem. Thus, sparsity of the difference to the prior image can be replaced by sparsity of the difference to the weighted mean:

χ t (f dyn ) = M -1 m=0 f m -f P 1 . ( 5
χ t (f dyn ) = M -1 m=0 f m - M -1 m ′ =0 w m ′ f m ′ 1 .
(5.14)

In the following, let us denote CS-sbs+t the reconstruction algorithm in which we use the temporal sparsity given by of Eq. (5.14) jointly with the SBS, and CS-wm+tv the reconstruction algorithm in which we use the same temporal sparsity jointly with TV minimization. All constraints must be satisfied at once. In absence of a direct expression for the proximal operator that is associated to CS reconstructions with multiple penalties, we turn to the iterative Dykstra-like proximal algorithm that was presented in Sec. 4.2.3 using:

• prox χ sbs that is the SBS operator of threshold λ sbs or prox χtv that is TV of weight λ tv and is handled by proximal splitting method such as the fast gradient projection [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF];

• and prox χt . Its direct expression is obtained by rewriting the temporal constraint as:

χ t = H(f dyn -f P ) 1 ,
where H is an orthonormal matrix. For the difference to the weighted mean, f P = 0 while H projects the temporal intensity at a given voxel into an average component and differential components that are thresholded. In the case of the prior image, H is set equal to the identity matrix.

In all our approaches, spatial hyperparameters λ sbs and λ tv are used with homotopy, while hyperparameters λ t can be assigned a fixed value without introducing a bias because temporal sparsity can be assumed.

Figure 5.9: Illustration of the proximal splitting for delayed opacification at stage s, that is the minimization of Eq. (5.12) using temporal sparsity of Eq. (5.14).

Weighted mean vs. prior image

We evaluated our approach on both simulated and clinical data.

Simulation study

In order to establish the performance of the CS reconstructions we simulated angiography data affected with motion by generating a phantom that consists of several phases of the same object with moving or varying intensity structures that are realistic with respect to the clinical data. To this aim, we built a numerical 2D phantom with three phases by adding simulated injected arteries (from 1500 to 6000 HU) to a 512 × 512 abdominal CT cross-section where background structures are valued between 1000 and 2000 HU. We simulated the acquisition in parallel geometry delivering 150 views in total of data with:

• a pulsatile motion of a few voxels only, which was generated using a real ECG signal;

• and a motion where each phase is defined by a subset of contiguous projections, resulting in a limited range sampling.

Projection simulation was carried out so as to mimic a dynamic acquisition: in practice, full scans are computed for each phase and a subset of views is selected along a given temporal signal.

Parameter Value 5.12 gives comparative assessment of the reconstruction quality. CS-sbs reconstruction is shown in Fig. 5.12.(a) to restore the temporal resolution of the vessels without introducing subsampling streaks, but the background is poorly depicted due to the substantial subsampling of each phase. Comparing to Fig. 5.12.(b), the prior image constraint of CS-piccs improved the background, while preserving the temporal resolution. A strong streak pattern remains in all RoIs however, that is a slightly attenuated version of the one present in the prior image (see Fig. Dealing with the limited range subsampling case whose results are displayed in Fig. 5.13, we observe that CS-sbs+t presented in Fig. 5.10.(c) outperformed CS-sbs and CS-piccs in the same manner. This visual inspection is confirmed by the final RMSD measurements1 d r given in Tab. 5.2 that are divided by two for the pulsatile case, and are comparable for the successive case. We quantified streak intensity by measuring the maximum intensity which is not a vessel in the RoI with dynamic injected vessels: for the pulsatile case, streak intensity decreased from 1872 HU for the prior image to 1636 HU and 1302 HU for CS-piccs and CS-sbs+t respectively; for the limited-angle case it decreased from 1954 HU to 1629 HU and 1406 HU. RMSD reference is taken as the reconstruction when no motion occurs during the scan acquisition.

Clinical Data

Results on real data are analyzed in Fig. 5.15 and Fig. 5.17 where standard clinical (i.e. static) reconstruction is compared to CS-piccs. The first dataset is the exam of the renal arteries (cf. Fig. 3.25), in which small displacements of the catheter (intense structure in the axial slices in Fig. 5.15) occur due to the pulsatile blood flow. We manually splitted the scan into four phases, relying on the vertical translation of the tip, that is visible on the 2D projections and plotted in Fig. left vessel in the axial slices in Fig. 5.17) did not appear opacified, while the left vertebral artery (LVA, most right vessel) was seen fully opacified during the whole scan. We thus splitted the scan into two phases that contained 72 and 75 projections respectively. Static reconstructions integrate all temporal variations: the position of the tip is blurred in the first exam (Fig. 5.15.(a)), while the intensity of the RVA peak is averaged to 9130 HU (Fig. 5.17.(a)), that is less than 60% of the LVA peak (15310 HU).

In addition, intense streak artifacts degrade the background: the displacements of the catheter yield a rotating pattern of positive and negative streaks in the first case, while the lack of opacification in the lateral projections yields horizontal and vertical negative streaks in the second case. CS-sbs+t reconstruction recovered some temporal resolution since temporal variations are visible in the reconstructed phases, even though our phase selection was approximate. This is particularly striking when looking at the associated maximum intensity projection (MIP) given in the second row of Fig. 5.17: the vertical translation of the catheter is well visible (Fig. 5.15.(b) and (c)); the RVA is only visible during phase 2, while the LVA is visible in both phases (Fig. 5.17.(b) and (c)). In terms of quantification the RVA peak was measured to be 1960 HU for phase 1 and 12290 HU for phase 2, while the LVA peak was measured to be around 13600 HU for both phases. As for the background CS-sbs+t reduced the streaks located close to the catheter tip (Fig. 5.15.(b) and (c)), and nearly eliminated the horizontal and vertical streaks of the RVA (Fig. 5.17.(b) and (c)) as confirmed by intensity profiles through the axial slice, that are plotted in Fig. 5.16 for the pulsatile motion case and in Fig. 5.18 for the limited-angle case.

Temporal Merge The sparsity of the temporal support is applied through a soft-thresholding operator. This therefore introduces a bias in the coupling between motion phases. Hence, λ t tuning results in a trade-off between background merge and temporal resolution. This trade-off is illustrated in Fig. 5.19: with λ t = 50 HU, that is displayed in Fig. 5.19(b), the artifact due to the lack of opacification of RVA is almost completely removed, but merging of the background is not satisfying compared to the back- ground given by FDK reconstruction shown in Fig. 5.19(a), while with λ t = 500 HU, that is displayed in Fig. 5.19(c), the artifact is more visible, but at the same time the background depiction is improved and very close to that of FDK.

Impact of the spatial penalty

In this section, we consider the same numerical 2D phantom that was used in Sec. 4.5.1. Five different phases of the phantom were produced by applying a global translation of 1 to 2 pixels along the x-axis and y-axis to the vessels. A detail of the simulated injected vessels on the slice background is given in Fig. 5.20 for three different motion phases. We simulated the acquisition in parallel geometry delivering 150 views in total of data with a pulsatile motion that results in a uniform subsampling for each motion phase.

Parameter Value τ 0.9 S 25 backgrounds I 1 iteration / background Table 5.3: Reconstruction parameters used for evaulating the impact of the spatial penalty. Figure 5.21 shows static FDK reconstruction as reference for the background quality. Note that despite the very small displacements of the vessels, the background of our simulation was severely de-teriorated. In particular the depiction of the bony structures is pretty poor. Furthermore, because many structures are not static, the blurred area is important. Figure 5.22 allows for comparative assessment of the reconstruction quality with no spatial constraints, SBS and TV. Each subfigure displays the reconstructed slice (Fig. 5.22(a)) and two regions of interest (RoIs) of the phantom with a modified windowing so as to enhance different structures of the image : the dynamic injected vessels over static background for two phases (Fig. 5.22(b)) and the static soft tissues (Fig. 5.22(c)). The introduction of the temporal sparsity assumption in the reconstruction, that is shown in the first row of Fig. 5.22, results in restoration of vessel temporal resolution accompanied by effective artifact reduction (bone depiction). However, significant streak artifacts remain, that prevent full restoration of the soft tissues. Excellent reconstruction results were obtained for CS-sbs+t and CS-wm+tv that are displayed in the second and third row of Fig. 5.22 respectively. Both effectively eliminated the motion artifacts, providing satisfying temporal resolution and background depiction. Note that the air-tissue interface that was blurred with static CS-sbs is now finely recovered with CS-sbs+t. It is difficult to conclude on the best approach between the SBS operator and TV based on the visual inspection of these reconstructions only. That is why we monitored the convergence of RMSD values2 for both CS-sbs+t and CS-wm+tv. Resulting log d r curves with respect to the number of iterations over the whole image are plotted in Fig.

CS reconstruction with motion modelling

The weighted mean is computed from different motion phases, thus contains motion artifacts and blur.

No merging is achieved in these areas. It is possible to go beyond simple sparsity and reduce the number of unknowns of the linear system given in Eq. (5.1) by incorporating a motion model that allows for describing the correlation between two phases, and thus reconstructing one single volume. The object of this contribution is to strengthen the proposed temporal constraint of CS-wm by computing a mean image that accounts for some (i.e. even incomplete) inter-phase motion compensation based upon an a priori motion knowledge. This is achieved by applying a transformation to deform the reconstruction of each phase f m to obtain a match with the reconstruction of arbitrarily chosen reference phase f r . Sparsity of the difference to the motion compensated weighted mean incorporates the motion information and can still be applied through a proximal operator. It follows CS A+T approach in which the following functional is minimized: fr,χ A +χ T λ A ,λ T = arg min

fr∈R K Q dyn (T f r ) + λ A • χ A (T f r ) + λ T • f r - M -1 m=0 w m T m f r 1 = χ T (f r )
.

(5.15)

We can distinguish two cases:

• the motion model is exact and complete: there is no guarantee that the data are fully sampled, so CS reconstruction can be applied to correct for remaining subsampling issues in the data;

• the motion model is locally exact but incomplete: the motion model is introduced to improve the CS constraints.

The ℓ 1 -norm handles the fact that locally the motion model might not be exact or complete and will not be used.

Non-rigid registration

The purpose of registration is to find a matching between phase f r and phase f m . The transformation T m expressed in Sec. 5.2 can be defined as the identity transformation I of which a dense displacement field D m has been added:

T m = I + D m .
The estimation of the optimal dense displacement field D m is carried out through non-rigid registration.

It is calculated through the optimization of a similarity criterion S that is used to estimate the quality of the registration: arg min Ψ S (T m f r -f m ), (5.16) where Ψ refers to the optimization strategy that allows for maximizing S.

Deformable registration methods can be generally categorized into two distinct groups:

• feature-based methods [START_REF] Davatzikos | Image registration based on boundary mapping[END_REF] that often rely on a preprocessing step to extract a set of salient features (e.g. points, edges) and establish correspondences on this basis. Since only the extracted features are used for registration, such methods are computationally very efficient. However, their performance is limited by the feature extraction;

• voxel-wise methods [START_REF] Zikic | Linear intensity-based image registration by markov random fields and discrete optimization[END_REF] that consists in maximizing a similarity criterion S based on voxel-wise attributes (e.g. intensities, distribution of the intensities).

In order to reduce the dimensionality of the problem, a current strategy consists in introducing a sparse set of control points and a set of interpolation functions that recovers the displacement specific to each voxel. Because the displacement of a control point has only local influence on the displacement field, cubic B-splines are often used as interpolation functions. The number of control points has to be a good compromise between speed and registration accuracy. That is the reason why a coarse-to-fine approach is often used in practice.

In this work, we investigated the use on angiographic data of two registration algorithms of which the basic principles are briefly presented hereafter. We rely on these two algorithms to highlight two key aspects of a registration algorithm: discriminative attributes and an efficient optimization.

DRAMMS Ou et al. [START_REF] Ou | Dramms: Deformable registration via attribute matching and mutualsaliency weighting[END_REF] developed the Deformable Registration via Attribute Matching and Mutual-Saliency weighting (DRAMMS) algorithm: it is an hybrid approach that alleviates the main drawbacks of both voxel-wise and feature-based approaches, and demonstrates large applicability with high robustness and accuracy, at the expense of a heavy computation. At each control point, a rich set of attributes (that was chosen as the multi-scale and multi-orientation Gabor attributes) is extracted then optimally selected, that allows for establishing strong correspondences between registered images. In addition, DRAMMS algorithm includes a weighting function called 'mutual-saliency' that reduces the impact of outlier regions by assigning weights that reflect the matching uniqueness between a pairs of voxels.

DROP Glocker et al. [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF][START_REF] Glocker | Deformable medical image registration: Setting the state of the art with discrete methods[END_REF] proposed the Deformable Registration using Discrete Optimization (DROP) algorithm: it is a computationally attractive intensity-based deformable registration that consists in recasting the registration problem as a Markov Random Field (MRF). In this context, the equivalent MRF problem accounts for energy minimization in a graph where nodes correspond to control points and a set of labels is associated with a set of deformations. Each control point is assigned a label such that the similarity measure between the source and the target is finally optimal for all voxels. Note that the optimization procedure is independent from the graph construction, and consequently any similarity criterion can be used. In this manuscript, because vascular structures cannot be recovered on the only basis of their intensities, we consider the sum of absolute differences plus the sum of gradient inner products (SADG) as a suitable similarity criterion:

S (T m f r -f m ) = 1 -γ K k |(f m ) k -(T m f r ) k )| + γ K k (∇f m ) k |(∇f m ) k | • (∇T m f r ) k ) |(∇T m f r ) k )| ,
where γ controls the balance between the sum of absolute differences and the sum of gradient inner products.

Performance of CS reconstructions using a prior deformation field

Reconstructions using a prior deformation field were assessed using a numerical 2D phantom with four phases by adding simulated injected arteries (from 1500 to 6000 HU) to a 512 × 512 abdominal CT crosssection where background structures are valued between 1000 and 2000 HU. We simulated the acquisition in parallel geometry delivering 150 views in total of data with a cyclic uniform motion. Reconstruction settings are given in Tab. 5.4. Here, we chose to rely on SBS to reduce angular streak artifacts, while estimating the prior deformation field by elastic registration of the phantom phases using DRAMMS algorithm3 . Figure 5.24 shows the deformation field (Fig. 5.24(a)) that was obtained when registering phase 1 (Fig. 5.24(b)) and phase 3 (Fig. 5.24(c)). Figure 5.25 shows reconstruction results using the prior deformation field. Application of the deformation field as post-processing to the reconstruction (i.e. elastic deformation of each phase to match the reference time-point, then computation of a weighted mean) recovers the vessel temporal resolution, but fails at removing streak artifacts that are associated to the subsampling of the moving structures. Hence the background of FBP reconstruction in identical. In the particular case of moving soft tissues (low contrasts), for which artifacts are more prone to be local, the approach may generate satisfying improvements, but such improvements no longer hold in the case of moving injected vessels (high contrasts). This can be perhaps more intuitively understood when assimilating streak artifacts to shadows in a photograph. The application of a deformation field that models the displacement of an object does not affect its shadow. The position of the light source and mechanism that leads the shadow formation are needed to generate a shadow that corresponds to the new position of the object. The same reasoning is valid in tomography. Fig. 5.25(c) displays CS A+T reconstruction result, that allows for mitigation of the streak artifacts. Computation of the RMSD over the whole slice confirms the significant image quality improvement that is provided by CS A+T approach: d r = 0.11 is found, which corresponds to about half the value that are found for static FBP reconstruction (d r = 0.21) and for FBP reconstruction on a single phase (d r = 0.24) that is not shown here.

Let us now consider the same renal exam that was used for producing the reconstruction results previously presented in Fig. 5.15. Since the registration task proved not being robust for processing subsampled volumes, in particular because of the prominent streak artifacts, the deformation field used as prior information in CS A+T reconstruction was calculated by registering previously obtained CS-Psbs reconstructions. . The motion of the catheter is satisfyingly recovered. However, the interpolation step applies the deformation to each phase introduced some smoothing of the reconstruction. Because the estimation of the computed deformation field contains many errors, this smoothing tends to spread the information within all volumes and slices displayed in Fig. 5.27 suffer from inter-slice diffusion (i.e. diffusion along the z-axis) of the rachis.

Discussion

We proposed a dynamic reconstruction that relies on a CS approach with spatio-temporal constraints.

The evaluation of the algorithm with numerical experiments and two typical angiographic datasets demonstrated qualitative and quantitative improvements. Limited assumptions were considered:

• high-intensity sparse structures but non-sparse background;

• phase temporal correlation. In order to simplify the presentation, the temporal constraint was also set equal for all voxels. This is not a requirement: in the case of a delayed opacification, each voxel could be treated independently, with its own phase selection and associated temporal constraint and proximal operator; • possible modelling of the vessel displacements through the estimation of the inter-phase deformation field. With the current implementation, such modelling is however not perfect for the still background. Future works shall aim at iteratively estimating the inter-phase deformation field. The computation at each iteration of the reconstruction process is likely to improve the quality of the reconstruction by refinement. One of the limitations to the introduction of a motion model into the reconstruction is however that variations in contrast-enhancement cannot be handled by elastic deformations, and hence, leads errors in the estimation of the deformation field. An adapted constraint relaxation approach could possibly compensate for the residual motions that are not modelled;

• similar attributes in the source and the target phases. In rotational angiography, phase-wise reconstructions account for the angular sampling and hence, artifacts may strongly differ from one phase to another. Thus, in the particular case of uniformly subsampled phases, the pattern of streak artifacts that equally deteriorate each phase is rotated from one phase to another. Because the streak artifacts deteriorate the entire image and have high intensities, it is generally not possible to perform a robust registration while handling them as outliers. When the phase sampling is very poor, it is even very likely that the anatomical structures result in less salient attributes than streaks themselves. Consequently, the estimation of the dense deformation field is a challenging task with state-of-the-art registration algorithms and a bottle-neck for reconstructions involving motion modelling. Possible improvements of the registration task include registering the moving (vascular) structures only [START_REF] Heibel | Discrete tracking of parametrized curves[END_REF] instead of registering the entire volume. This would require the proper segmentation of the contrast-injected vessels. To this aim, it would be possible to use CS reconstruction that we discussed in the previous chapter since it enables a thresholded segmentation of the sparse structures and then decrease the threshold level by continuation. Another possible strategy could consist in using 'simplified images' that are free of subsampling streak artifacts instead of the original ones. Such images could be obtained through a TV-penalized reconstruction for instance;

• prior phase selection. Reconstructions were obtained in a clinical context with manual phase selection. For an actual clinical usage, automatic phase selection is a requirement. It is challenging, but to a lesser degree than motion estimation and modelling.

Chapter 6

Conclusion 6.1 Outcome

This work tackled the reconstruction of tomographic acquisitions that were acquired with C-arm systems in the context of interventional radiology. It deals more specifically with the correction of motion artifacts, that originate from the temporal variations of the contrast-enhanced vessels and may severely deteriorate the overall image quality, and highlights a central aspect of tomography, that is data (angular) sampling.

The sampling trade-off inherent to the interventional context was discussed and shown to depend on both the C-arm technology, that was originally designed for real-time guidance, and on the angiographic procedures. Thus, the choice of the rotation speed depends on the clinical protocol: the rotation of the Carm gantry is used at maximal speed for injected acquisitions, while kept slow for non-angiographic exam in order to maximize the sampling. Sampling artifacts specifically induced by C-arm systems involve: the cone-beam geometry, that results in cone-beam artifacts, and the relatively low detector framerate, that results in uniform subsampling. This limits the use of three-dimensional imaging in interventional procedures. On the other hand, better soft tissue depiction would benefit from the high spatial resolution of C-arm systems, and thus participate to the development of novel clinical applications. The developments of the compressed sensing theory have shown that least-square minimization combined with sparse constraints expressed as ℓ 1 -penalties may yield satisfying data recovery in cases where sampling is not favourable, under well-defined assumptions. Although the compressed sensing assumptions are not strictly valid for X-ray tomography, its first application -exact reconstruction of a piecewise constant phantom from poorly angularly sampled data that simulated the X-ray sampling processopened tremendous perspectives for image quality improvement in X-ray tomography. Reconstructions of dynamic data based on compressed sensing were proposed in the literature. A particularly seducing CS approach is the PICCS algorithm that relies on using a sparse constraint based on a prior image -that may be the static filtered backprojection reconstruction-, so that the reconstruction, schematically, determines motionless areas from the complete acquisition, while motion-blurred areas are determined from reconstruction of well-identified subsampled motion phases only. We however observed that the PICCS algorithm is not applicable to intraarterial contrast injection, since intense streaks due to the presence of inconsistent contrast-enhanced vessel projections in the data makes the static reconstruction a poor prior image. In this line, the present work proposed an application of compressed sensing and a generalization of the PICCS algorithm that are adapted to rotational angiography with intraarterial injection.

For easily adding a series of ℓ 1 or ℓ 2 constraints to iterative least-square the proximal splitting framework that relies on proximal operators was used. An iterative implementation based on filtered backprojection and proximal splitting was proposed, that allows for computational efficiency. Several variants of iterative filtered backprojection were derived to deal with sparse penalties:

• soft background subtraction that ensures sparsity of the image ℓ 1 -norm;

• total variation that ensures sparsity of the gradient ℓ 1 -norm;

• non linear ℓ 2 -minimization such as anisotropic filtering. This work also emphasizes the fact that compressed sensing reconstruction is not reduced to ℓ 1 -norm minimization: starting from ℓ 1 -norm minimization and using a homotopy continuation strategy, relation-CHAPTER 6. CONCLUSION ship with empirical orthogonal matching pursuit approaches are observed, leading to the key contribution of this thesis that is the importance of homotopy with respect to regularization.

Compressed sensing reconstruction was applied to well-identified clinical cases with intraarterial injections and resulting sampling issues, where image quality improvement would be welcome. This work focused more specifically on improvement for both cerebral and abdominal imaging. In a static interventional context, artifact reduction was shown for clinical data in spite of data truncation in the following identified cases:

• cone-beam artifact correction;

• subtracted reconstruction;

• contrast-enhanced reconstruction;

• limited-angle reconstruction, but with a limited impact because the non-sparse background can only be reconstructed from the entire scan.

However, streak artifacts originated from metal or objects outside the FoV could not be removed. Further exploring compressed sensing for the dynamic interventional context, the following improvements were proposed:

• an original compressed acquisition design for digital subtracted rotational angiography;

• the introduction of an additional temporal sparse penalty within a multiple penalties framework for combining the compressed sensing reconstruction of sparse structures with the recovery of a static background with no motion assumption.

Perspectives

Several directions for future works were identified, that would either be required for industrialization of this work or extend its clinical applications:

• For an actual clinical usage of the dynamic approach, automatic phase selection from the whole acquisition is a requirement. It is challenging, but to a lesser degree than motion estimation and modelling;

• Still for the dynamic approach, the stagewise orthogonal matching can be derived in a temporal and iterative variant as an alternative to ℓ 1 -penalization. In addition, the constraint on the temporal sparsity allows for obtaining partial motion knowledge and consequently for revisiting motion compensated approaches, since the exact and complete knowledge of the motion is no longer required. The introduction of a motion model was briefly discussed and shown to have potential to either complement (considering a cardiac motion) or supplement (considering a respiratory motion) the compressed sensing reconstruction. Further investigations on clinical data are required to evaluate these preliminary results;

• Future work shall aim at applying the proposed mathematical framework to a larger class of algorithms and thus associate iterative filtered backprojection with standard image processing tools that cannot be described through a variational formulation (e.g. bilateral filtering, diffusion filtering, NL-means);

• This work has potential application to CT (e.g. helical artifact correction, dynamic reconstruction of cardiac and perfusion data). This implicitly raises the need for comparing the compressed sensing approach (low angular sampling strategy) with a noise modelling approach that is coupled with a sparse temporal regularization (dose reduction strategy). As far as noise can be modelled through a weighted least-square, the reconstruction framework presented here is still valid. Compressed sensing reconstruction is however no longer suitable for dealing with Poisson models and consequently alternative approaches must be considered [START_REF] Pustelnik | Méthodes proximales pour la résolution de problèmes inverses: application à la tomographie par émission de positrons[END_REF][START_REF] Ting | Sparse Image Reconstruction for Molecular Imaging[END_REF].

  2.19 Helical CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 C-arm systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 C-arm systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Digital flat-panel detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Example of rotational scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Coordinate systems in cone-beam geometry . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Cone-beam geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Coordinate systems in fan-beam geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Scatter reduction through cone angle reduction . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Truncation correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Comparative evaluation of the image quality between C-arm systems and CT scanners . . 3.11 Overview of the standard interventional workflow . . . . . . . . . . . . . . . . . . . . . . . 3.12 Fluoroscopic 2D views for digital subtraction angiography . . . . . . . . . . . . . . . . . . 3.13 Digital subtraction rotational angiography reconstruction . . . . . . . . . . . . . . . . . . 3.14 Example of aneurysm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Coil embolisation of an aneurysm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Roadmapping by merging 3D with fluoroscopy . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Three-dimensional imaging for selective chemoembolisation of a liver tumour . . . . . . . 3.18 Vascular anatomy of the hepatic artery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Follow-up CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Example of subsampling artifacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Effect of respiration on the performance of standard FDK reconstruction . . . . . . . . . . 3.22 Example of artifacts due to small displacements of contrast-enhanced vessels . . . . . . . 3.23 Delayed opacification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 Example of artifacts caused by contrast enhancement variations . . . . . . . . . . . . . . . 3.25 Angiography data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 Evaluation of the dynamic of the motion in a scan . . . . . . . . . . . . . . . . . . . . . .

Figure 2 .

 2 Figure 2.1: Coordinate systems in parallel-beam geometry.

Figure 2 . 2 :

 22 Figure 2.2: (a) Imaged object and (b) its sinogram.

Figure 2 . 3 :

 23 Figure 2.3: Relationship between spatial, Fourier and Radon domains.

0 +∞ 0 [F

 00 k y )e i2π(kxx+kyy) dk x dk y By applying the change of variables from Cartesian coordinates (k x , k y ) into polar coordinates (k u , φ), it follows that:f (x, y) = 2π 2D f ](k u cos φ, k u sin φ)e i2πku(x cos φ+y sin φ) k u dk u dφ

[F[F

  1D Rf ](k u , φ)e i2πku(x cos φ+y sin φ) k u dk u dφ= 1D Rf ](k u , φ)|k u |e i2πkuu dk u dφ = π 0 [p P B,φ * h](u)dφ ,

Figure 2 . 4 :

 24 Figure 2.4: Band-limited ramp filter: (a) frequency response, (b) impulse response.

Figure 2 . 5 :

 25 Figure 2.5: Angular sampling patterns. (a) Full sampling. (b) Uniform subsampling.

Figure 2 . 6 :

 26 Figure 2.6: Impact of the number of integral lines and bins on FBP reconstruction. (a) Imaged object. (b) FBP reconstruction from 64 integral lines that are equiangularly spaced over 180 • and 256 bins in parallel geometry. (c) FBP reconstruction from 256 integral lines that are equiangularly spaced over 180 • and 64 bins in parallel geometry.

Figure 2 . 7 :

 27 Figure 2.7: Angular subsampling and streak artifacts. (a) FBP reconstruction from 32 integral lines that are equiangularly spaced over 180 • in parallel geometry. (b) Difference between image (a) and the image object previously displayed in Fig. 2.6(a).

Figure 2 . 8 :

 28 Figure 2.8: (a) Limited angle sampling subset. (b) Nonuniform subsampling subset.

Figure 2 . 9 :

 29 Figure 2.9: Impact of angular coverage on FBP reconstruction. (a) Phase 1 and (b) phase 2 of an imaged object with background of value 1000 and moving structures of contrast 5000 -value range: 0 to 3000. Simulation of a parallel-beam sinogram with 300 equiangularly spaced integral lines over 180 • such that the imaged object is seen in phase 1 over the 200 first angles, while in phase 2 over the last 100 angles. Partial FBP reconstructions: (c) phase 1 and (d) phase 2. (e) static FBP reconstruction.

Figure 2 . 10 :

 210 Figure 2.10: Motion artifacts with low contrast structures (a) Imaged object with background of value 1000 and moving structures of contrast 100 -display range: 900 to 1110. Simulation of a parallel-beam sinogram with 300 equiangularly spaced integral lines over 180 • . Uniform periodic displacement of the moving structure: (b) sinogram; (c) static FBP reconstruction.

Figure 2 . 11 :

 211 Figure 2.11: Motion artifacts with high contrast structures and dependency with respect to sampling. (a) Imaged object with background of value 1000 and moving structures of contrast 5000 -display range: 0 to 2000. Simulation of a parallel-beam sinogram with 300 equiangularly spaced integral lines over 180 • . Uniform periodic displacement of the moving structure: (b) sinogram; (e) static FBP reconstruction. Nonuniform periodic displacement of the moving structure: (c) sinogram; (f) static FBP reconstruction. Sudden intensity change: (d) sinogram (g) static FBP reconstruction.

Figure 2 .

 2 Figure 2.12: (a) Schematic representation of the elements of an X-ray tube [1]. (b) Photograph of an X-ray tube (GE Healthcare, Chalfont St. Giles, UK).

Figure 2 . 14 :

 214 Figure 2.14: Photograph of a GE Lightspeed (GE Healthcare, Chalfont St. Giles, UK).

Figure 2 . 15 :

 215 Figure 2.15: Radiologic terminology and conventions. Axial planes run from the cranial end to the caudal end. Coronal planes run from the ventral side to the dorsal side. Sagittal planes runs from the left side to the right side or vice-versa. In radiology, by convention structures are observed from the caudal end up with the patient in front of the observer and hence, left structures of the body are represented on the right side in axial planes (and vice-versa).

Figure 2 . 16 :Figure 2 . 17 :

 216217 Figure 2.16: Illustration of the windowing principle to modify the dynamic of the displayed image. (a) A wide window is chosen to display tissues largely differing in density: WL = 35 HU and WW = 2000 HU. (b) A narrow window is chosen to visualize small contrasts: WL = 35 HU and WW = 100 HU

  (2.19) shows that if β ∈ [0, 2π], all measurements γ will sample interval[0, 2π], while if β ∈ [0, π], detector bin γ = -Γ samples [-Γ, π -Γ], bin γ = 0 samples [0, π] and bin γ = Γ samples [Γ, π + Γ]. Thus sufficient condition for sampling [0, π] for all γ is that β ∈ [-Γ, π + Γ].It does not provide a one-to-one correspondence with parallel-beam measurements since some points are sampled twice while others are sampled only once. In practice the 2π coverage is preferred since it provides a uniform double sampling of the parallel-beam geometry.

Figure 2 . 18 :

 218 Figure 2.18: Short scan. Sinogram measurement plotted in the (β,γ) domain.

Figure 2 .

 2 Figure 2.19: Helical CT.

Figure 3

 3 Figure 3.1: C-arm systems: (a) mobile C-arm (GE Healthcare, Chalfont St. Giles, UK); (b) biplane system (GE Healthcare, Chalfont St. Giles, UK); (c) Artis zeego (Siemens Healthcare, Erlangen, Germany); (d) GE IGS 730 (GE Healthcare, Chalfont St. Giles, UK).

Figure 3 . 2 :

 32 Figure 3.2: Photograph of an Innova system (GE Healthcare, Chalfont St. Giles, UK). Elements of the imaging system and their articulations.

Figure 3 . 3 :

 33 Figure 3.3: (a) Schematic representation of the layers that compound the sensor used in flat-panel detectors. The detector is based on a crystalline CsI scintillator and an a-Si diode. Active readout matrix includes readout electronics based on a small thin-film transistor. (b) Photograph of the flat panel detector of Innova 4100 with its packaging. (c) Photograph of the flat panel detector of Innova 4100.

Figure 3 . 4 :

 34 Figure 3.4: Example of rotational scan. Each view corresponds to a particular position of the C-arm. Vertical arrow: number of slices. Horizontal arrow: axial field of view.

Figure 3 . 5 :

 35 Figure 3.5: Coordinate systems in cone-beam geometry.

Figure 3 .

 3 6 illustrates the different coordinate systems and necessary transforms to geometrically model a cone-beam projection. Let us denote (x ′ S , y ′ S , z ′ S ) the absolute position of source S, and respectively (x ′

Figure 3 . 7 :

 37 Figure 3.7: Coordinate systems in fan-beam geometry with a flat detector.

Figure 3 . 8 :

 38 Figure 3.8: (a) Standard reconstruction. (b) Scatter reduction through cone angle reduction.

Figure 3

 3 Figure 3.9: (a) Standard reconstruction. (b) Truncation correction.

Figure 3 . 10 :

 310 Figure 3.10: The comparative evaluation of three-dimensional reconstructions of neurological data from (a) a C-arm system and (b) a CT scanner showing axial slices illustrates superior image quality achieved with the CT scanner. In particular, differentiation of the brain soft tissues is easier due to higher contrast resolution. (c) Contrast resolution with a C-arm system.

Figure 3 . 11 :

 311 Figure 3.11: Overview of the standard interventional workflow.

Figure 3 . 12 :

 312 Figure 3.12: Fluoroscopic 2D views for digital subtraction angiography. (a) Mask 2D view. (b) Contrast 2D view. (c) Corresponding subtracted 2D view.

Figure 3 . 13 :

 313 Figure 3.13: Digital subtraction rotational angiography reconstruction. First row: transaxial slice. Second row: MIP visualization. (a) Mask volume. (b) Contrast volume. (c) Subtracted volume

Figure 3 . 14 :

 314 Figure 3.14: Example of saccular aneurysm (volume rendering).

Figure 3 . 15 :

 315 Figure 3.15: Coil embolisation of an aneurysm. (a) A catheter is advanced to the location of the aneurysm. (b) Detachable coils are inserted through the catheter. (c) Coils completely fill the aneurysm sac so that the blood flow can no longer be directed to the aneurysm.

Figure 3 . 16 :

 316 Figure 3.16: Roadmapping by merging 3D with fluoroscopy. (a) Fluoroscopy. (b) 3D vasculature. (c) Fluoroscopy merged with 3D.

Figure 3 . 18 :

 318 Figure 3.18: Vascular anatomy of the hepatic artery.

Figure 3 .Figure 3 . 20 :

 3320 Figure 3.19: (a) Innova 3D image performed during the procedure. (b) Embolisation results (indicated by orange arrows) are demonstrated by post-procedural CT image.

Figure 3 . 21 :

 321 Figure 3.21: Effect of respiration on the performance of standard FDK reconstruction. Transaxial slice of a swine abdominal anatomy. First row: FDK reconstruction from a free-motion scan. Second row: FDK reconstruction from a scan acquired during free breathing. (a) Transaxial view. (b) Coronal view. (c) Sagittal view.

Figure 3 . 22 :Figure 3 . 23 :

 322323 Figure 3.22: Example of artifacts due to small displacements of contrast-enhanced vessels. (a) Transaxial slice is deteriorated by motion streak artifacts. (b) MIP image is more robust to motion streaks.

Figure 3 . 24 :

 324 Figure 3.24: Example of artifacts caused by contrast enhancement variations. Delayed opacification: (a) transaxial slice, (b) coronal slice. Incomplete dilution: (c) transaxial slice, (d) coronal slice.

  .26. Note that the axial component of the catheter motion is somehow visible in the sinogram representations of the projection data: it results in oscillation in Fig. 3.25(b), Fig. 3.25(c) and Fig. 3.25(d).

Figure 3 . 25 :

 325 Figure 3.25: Angiography data. (a) 2D view. (b), (c) and (d) Sinogram views that correspond to different horizontal lines of (a) and where oscillations due to the variations in contrast can be seen.

Figure 3 . 26 :

 326 Figure 3.26: Evaluation of the dynamic of the motion in a scan. Manual identification of the temporal signal that is associated to the motion of the catheter.

Figure 4 . 1 :

 41 Figure 4.1: Geometric interpretation of the ℓ 1 and ℓ 2 norms. (a) ℓ 1 ball: the anisotropy of the ℓ 1 ball favors sparse solution. (b) ℓ 2 ball: minimizing the ℓ 2 norm generally does not allow for recovery of sparse solutions.

Figure 4 . 2 :

 42 Figure 4.2: Illustration of OMP principle. (a) Sparse phantom. (b) FDK reconstruction from 10 views. (c) Detected components after thresholding. (d) FDK reconstruction of the residual. (e) Detected components after thresholding. (f) Final reconstruction: sum of (c) and (e).

Figure 4 . 3 :

 43 Figure 4.3: Visual comparison of (top row) FDK and (bottom row) iFDK+ with regard to the sampling rate. WL=1000 HU -WW=2000 HU. (a) Full scan (about 150 views). (b) Subsampling factor: 2. (c) Subsampling factor: 4. (d) Subsampling factor: 8.

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Comparison of FDK and iFDK+ with regard to the sampling rate. Intensity profile curves drawn through the axial slice presented in Fig. 4.3 for (a) the full scan (about 150 views) (b) a subsampling factor equal to 8. Reconstruction using FDK is plotted in blue, while reconstruction using iFDK+ is plotted in red.

. 36 )

 36 while initializing the computation of fχ sbs λ (s) sbs , i.e. minimization of Eq. (4.36) at stage s with the solution computed at the previous stage fχ sbs λ (s-1) sbs . The minimization at λ sbs = λ (s)

. 37 )Figure 4 . 6 :

 3746 Figure 4.6: Sparsity in the gradient domain. (a) Gradient image of a C-arm CT slice. (b) Corresponding intensity histogram.

Figure 4 . 7 :

 47 Figure 4.7: Sparsity in the wavelet domain (Daubechies 9/7). (a) Reference image. (b) Compression via removal of 90% of the coefficients of the first detail subband. (c) Compression via removal of 90% of the coefficients of the two first detail subbands.

Figure 4 .Figure 4 . 8 :

 448 Figure 4.8: 2D numerical static phantom. (a) Full slice. HU range: 340 to 1900. (b) Detail of the injected vessels. HU range: 850 to 1230. (c) Detail of the soft-tissues. HU range: 600 to 1200.

Figure 4 . 9 :

 49 Figure 4.9: Reconstruction results with positivity and with regularization strategies. First row: reconstructed image. Second row: zoomed detail with HU range from 850 to 1230. Third row: zoomed detail with HU range from 600 to 1200. Reconstruction result of: (a) experiment n • 1. (b) experiment n • 2. (c) experiment n • 3. (d) experiment n • 4. (e) experiment n • 5.

3 :

 3 Generic parameters of iFBP-derived algorithms.

Figure 4 .

 4 Figure 4.11 compares CS reconstruction results relying on different spatial sparsity assumptions.Subsampling resulted in an iFBP+ reconstruction (Fig.4.11(a)) that was deteriorated by a rotating pattern of streaks, of which the strongest ones come from the injected vessels and the bony structures. Soft background subtraction (SBS) corrected for streak artifacts coming from highly intense structures,

Figure 4 . 10 :

 410 Figure 4.10: Reconstruction results with homotopy strategies. First row: reconstructed image. Second row: zoomed detail with HU range from 850 to 1230. Third row: zoomed detail with HU range from 600 to 1200. Reconstruction result of: (a) experiment n • 6. (b) experiment n • 7. (c) experiment n • 8. (d) experiment n • 9. (e) experiment n • 10.

w = 1000 to 0 Table 4 . 4 :

 044 Parameters of the sparsity constraint.

Figure 4 . 11 :

 411 Figure 4.11: Subsampled reconstructions. First row: reconstructed image. Second row: zoomed detail on injected vessels. Third row: zoomed detail on soft tissues. (a) iFBP+ reconstruction. (b) CS-sbs reconstruction. (c) CS-tv reconstruction. (d) CS-w reconstruction.

Figure 4 . 12 :

 412 Figure 4.12: Convergence curves computed on the reconstructions displayed in Fig. 4.11.

Figure 4 . 13 :

 413 Figure 4.13: iFBP with anisotropic diffusion. The strength of anisotropy is defined by the full width at half maximum (fwhm). First row: reconstructed image. Second row: zoomed detail on injected vessels. Third row: zoomed detail on soft tissues. (a) Reconstruction with fwhm = 10. (b) Reconstruction with fwhm = 1. (c) Reconstruction with fwhm = 10 to 1.

Figure 4 .Figure 4 . 14 :

 4414 Figure 4.14: Visual comparison between FDK, positivity constraint and SBS in the subtracted case. HU range: 550 to 1050. (a) FDK reconstruction. (b) iFDK+ reconstruction. (c) CS -sbs reconstruction.

Figure 4 . 15 :Figure 4 . 16 :Figure 4 . 17 :

 415416417 Figure 4.15: Intensity profile curves through the slice shown in Fig. 4.14.

Figure 4 . 18 :

 418 Figure 4.18: Intensity profile curves for comparison of positivity constraint and SBS in the non-subtracted case: (a) through the axial slices shown in Fig. 4.16, (b) through the sagittal slices shown in Fig. 4.17.

Figure 4 . 19 :Figure 4 . 20 :Figure 4 . 21 :Figure 4 . 22 :Figure 4 . 23 :

 419420421422423 Figure 4.19: CS reconstruction and metal artifacts. HU range: 750 to 1350. First row: iFDK+ reconstruction. Second row: CS -sbs reconstruction. Third row: CS -tv reconstruction. Sampling rate: (a) 600 views. (b) 150 views. (c) 75 views.

Figure 4 .Figure 4 . 24 :

 4424 Figure 4.24: Reconstruction results in the subtracted case (HU range: -1000 to 6000). (a) FDK reconstruction of 150 views over 0-200 • . Reconstruction of 75 views over 100-200 • : (b) fiFDK+ reconstruction, and (c) fCS-sbs reconstruction.

Figure 4 . 25 :

 425 Figure 4.25: Profile curve through the slice shown in Fig. 4.24.

Figure 4 . 26 :

 426 Figure 4.26: Non subtracted RVA (HU range: 0 to 2400). (a) FDK reconstruction of 150 views over 0-200 • . (b) FDK reconstruction of 75 views over 100-200 • . (c) CS-sbs reconstruction of 75 views over 100-200 • .

Figure 5 . 1 :

 51 Figure 5.1: Dynamic data acquisition model: the temporal variations encountered during the acquisition within a 2-phase object are modelled as a succession of phase f 0 (in light gray) and phase f 1 (in dark gray).

Figure 5 . 2 :

 52 Figure 5.2: Toy example. (a) 2D numerical phantom. (b) Simulation of the sinogram acquired in parallel geometry with 150 integral lines while the phantom is rotating at the same speed as the system gantry. (c) FBP reconstruction. (d) Motion-compensated reconstruction.

Figure 5 . 3 :

 53 Figure 5.3: DSRA acquisition protocol design. The mask angular sampling pattern is plotted with dashed gray arrows, while the contrast angular sampling pattern is plotted with dashed black arrows. (a) Homogeneous acquisition: contrast and mask scans sample the same set of equiangularly spaced positions. (b) Heterogeneous acquisition with δ = 0.5∆φ: contrast and mask scans sample two sets of equiangularly spaced interleaved positions.

Figure 5 . 4 :

 54 Figure 5.4: Impact of the acquisition protocol on the subtracted volume (N = 75 views). (a) FDK reconstruction for the homogeneous protocol: the redundant background structures and the associated streaks are removed. (b) FDK reconstruction for a heterogeneous protocol (δ = 0.5∆φ): streaks associated to the background structures cannot be removed. Vessel streaks remain with both protocols.

Figure 5 . 5 :

 55 Figure 5.5: Background sampling improvement with a heterogeneous protocol: each scan samples N = 75 views and angular shift δ = 0.5∆φ. (a) FDK reconstruction of the mask volume. (b) FDK reconstruction of the contrast volume. (c) A posteriori combination of the mask and the contrast volumes: the background sampling is doubled, but vessel streaks remain.

Figure 5 . 6 :

 56 Figure 5.6: DSRA reconstruction from simulated data. HU range: 1020 to 1100. (a) Mask iFBP+ reconstruction with double sampling (300 views). (b) Mask iFBP+ reconstruction. (c) Contrast iFBP+ reconstruction with simulated parenchyma and injected vessels. (d) Mask CS t reconstruction with λ t = 5 HU. (e) Mask CS sbs+t reconstruction with λ t = 5 HU.

Figure 5 . 7 :

 57 Figure 5.7: DSRA reconstruction from clinical data. First line: Mask iFDK+ reconstruction with spatial positivity constraint. Second line: Mask CS sbs+t reconstruction with λ t = 5 HU. (a) Axial slice (HU range: 250 to 5350). (b) Detail of the petrous part of the left temporal bone in (a). (c) Subtraction image for an axial slice higher in the brain (HU range: -500 to 1500).

Figure 5 . 8 :

 58 Figure 5.8: Sparsity in the temporal domain. (a) Zoomed detail of the MIP representation for two time points with small displacement of a sparse structure over a static background. The position of the object in the MIP corresponding to time point 1 (top image) is displayed in orange and reported into the MIP corresponding to time point 2 (bottom image). (b) Corresponding intensity histogram.

  λ t = 50 HU Table5.1: Reconstruction parameters used for comparing temporal penalties.

Figure 5 .

 5 Figure 5.10 displays two regions of interest (RoI) of the 2D numerical phantom: with static structures only (Fig. 5.10.(a)) and with simulated dynamic injected vessels over static background for two phases (Fig. 5.10.(b) and Fig. 5.10.(c)).Figure5.12 gives comparative assessment of the reconstruction quality. CS-sbs reconstruction is shown in Fig.5.12.(a) to restore the temporal resolution of the vessels without introducing subsampling streaks, but the background is poorly depicted due to the substantial subsampling of each phase. Comparing to Fig.5.12.(b), the prior image constraint of CS-piccs improved the background, while preserving the temporal resolution. A strong streak pattern remains in all RoIs however, that is a slightly attenuated version of the one present in the prior image (see Fig. 5.11.(a) and see Fig.5.11.(b)). The best restoration is thus obtained by CS-sbs+t as shown in Fig. 5.12.(c). Dealing with the limited range subsampling case whose results are displayed in Fig. 5.13, we observe that CS-sbs+t presented in Fig. 5.10.(c) outperformed CS-sbs and CS-piccs in the same manner. This

Figure

  Figure 5.10 displays two regions of interest (RoI) of the 2D numerical phantom: with static structures only (Fig. 5.10.(a)) and with simulated dynamic injected vessels over static background for two phases (Fig. 5.10.(b) and Fig. 5.10.(c)).Figure5.12 gives comparative assessment of the reconstruction quality. CS-sbs reconstruction is shown in Fig.5.12.(a) to restore the temporal resolution of the vessels without introducing subsampling streaks, but the background is poorly depicted due to the substantial subsampling of each phase. Comparing to Fig.5.12.(b), the prior image constraint of CS-piccs improved the background, while preserving the temporal resolution. A strong streak pattern remains in all RoIs however, that is a slightly attenuated version of the one present in the prior image (see Fig. 5.11.(a) and see Fig.5.11.(b)). The best restoration is thus obtained by CS-sbs+t as shown in Fig. 5.12.(c). Dealing with the limited range subsampling case whose results are displayed in Fig. 5.13, we observe that CS-sbs+t presented in Fig. 5.10.(c) outperformed CS-sbs and CS-piccs in the same manner. This

  Figure 5.10 displays two regions of interest (RoI) of the 2D numerical phantom: with static structures only (Fig. 5.10.(a)) and with simulated dynamic injected vessels over static background for two phases (Fig. 5.10.(b) and Fig. 5.10.(c)).Figure5.12 gives comparative assessment of the reconstruction quality. CS-sbs reconstruction is shown in Fig.5.12.(a) to restore the temporal resolution of the vessels without introducing subsampling streaks, but the background is poorly depicted due to the substantial subsampling of each phase. Comparing to Fig.5.12.(b), the prior image constraint of CS-piccs improved the background, while preserving the temporal resolution. A strong streak pattern remains in all RoIs however, that is a slightly attenuated version of the one present in the prior image (see Fig. 5.11.(a) and see Fig.5.11.(b)). The best restoration is thus obtained by CS-sbs+t as shown in Fig. 5.12.(c). Dealing with the limited range subsampling case whose results are displayed in Fig. 5.13, we observe that CS-sbs+t presented in Fig. 5.10.(c) outperformed CS-sbs and CS-piccs in the same manner. This

  Figure 5.10 displays two regions of interest (RoI) of the 2D numerical phantom: with static structures only (Fig. 5.10.(a)) and with simulated dynamic injected vessels over static background for two phases (Fig. 5.10.(b) and Fig. 5.10.(c)).Figure5.12 gives comparative assessment of the reconstruction quality. CS-sbs reconstruction is shown in Fig.5.12.(a) to restore the temporal resolution of the vessels without introducing subsampling streaks, but the background is poorly depicted due to the substantial subsampling of each phase. Comparing to Fig.5.12.(b), the prior image constraint of CS-piccs improved the background, while preserving the temporal resolution. A strong streak pattern remains in all RoIs however, that is a slightly attenuated version of the one present in the prior image (see Fig. 5.11.(a) and see Fig.5.11.(b)). The best restoration is thus obtained by CS-sbs+t as shown in Fig. 5.12.(c). Dealing with the limited range subsampling case whose results are displayed in Fig. 5.13, we observe that CS-sbs+t presented in Fig. 5.10.(c) outperformed CS-sbs and CS-piccs in the same manner. This

Figure 5 . 10 :

 510 Figure 5.10: Zoomed details of the phantom used for comparing temporal penalties. (a) RoI 1 that contains static structures -HU display range: 800 to 1400. RoI 2 that contains dynamic injected vessels over static background for (b) phase 1 and (c) phase 2 respectively -HU display range: 100 to 2800.

Figure 5 .

 5 Figure 5.14 contains the plots of log d r with respect to the number of iterations for two RoIs: one covering the area with moving vessels displayed in Fig. 5.10.(b) and shown in graph 5.14.(a), the other covering the background displayed in Fig. 5.10.(a) and shown in graph 5.14.(b). Discrepancy between the prior image and the reference image over each RoI is shown as a black dashed line. The final deviations of graph 5.14.(a) reflect the improved recovery of the moving structures with all CS dynamic reconstructions. In graph 5.14.(b), CS-sbs high deviation reflects poor background recovery. For CSpiccs, background recovery nears the level of the prior image, but cannot improve it by construction. On the contrary, CS-sbs+t restores the background below the discrepancy level of the prior image.

3 . 26 .Figure 5 . 11 :

 326511 Figure 5.11: Static reconstructions. Pulsatile case: (a) RoI 1 that contains static structures -HU display range: 800 to 1400; (b) RoI 2 that contains dynamic injected vessels over static background for phase 1 -HU display range: 100 to 2800. Limited-angle case: (c) RoI 1 that contains static structures; (d) RoI 2 that contains dynamic injected vessels over static background for phase 1.

Figure 5 . 12 :Figure 5 . 13 :Figure 5 . 14 :

 512513514 Figure 5.12: Comparison of temporal penalties for the pulsatile case. First row: RoI 1 that contains static structures -HU display range: 800 to 1400. Second and third rows: RoI 2 that contains dynamic injected vessels over static background for phase 1 and phase 2 respectively -HU display range: 100 to 2800. Reconstruction algorithm: (a) CS-sbs; (b) CS-piccs; (c) CS-sbs+t.

Figure 5 . 15 :

 515 Figure 5.15: Reconstruction results for the pulsatile motion case. Axial slice and MIP representation. First row: axial slice detail -HU display range: 350 to 1550. Second row: MIP detail -HU display range: 1100 to 4400. Algorithms: (a) static reconstruction. Orange arrow indicates streaks that originate from the catheter tip; CS-sbs+t for (b) phase 1 and (c) phase 2. (d) and (e) Different parameters.

Figure 5 . 16 :

 516 Figure 5.16: Reconstruction results for the pulsatile motion case. Intensity profile curves drawn through the axial slices shown in Fig. 5.15: static reconstruction (black dashed line), CS-sbs+t for phase 1 (blue plain line) and phase 2 (red plain line).

Figure 5 . 17 :

 517 Figure 5.17: Clinical data limited-angle case. First row: axial slice detail -HU display range: (1) 350 to 1550 (2) 0 to 2400. Second row: MIP detail -HU display range: (1) 1100 to 4400 (2) 1000 to 9500. Algorithms: (a) static reconstruction. Orange arrows indicate horizontal and vertical streaks; CS-sbs+t for (b) phase 1 and (c) phase 2.

Figure 5 . 18 :

 518 Figure 5.18: Reconstruction results for the late opacification case. Intensity profile curves drawn through the axial slices presented in Fig. 5.17: static reconstruction (black dashed line), CS-sbs+t for phase 1 (blue plain line) and phase 2 (red plain line).

Figure 5 . 19 :

 519 Figure 5.19: Trade-off between background merge and temporal resolution. (a) Static reconstruction. (b) λ t = 50 HU. (c) λ t = 500 HU.

Figure 5 . 20 :

 520 Figure 5.20: 2D dynamic numerical phantom. Detail of the injected vessels for (a) phase 1, (b) phase 2 and (c) phase 3. HU range: 850 to 1230.

Figure 5 . 21 :Figure 5 . 22 :

 521522 Figure 5.21: Static FBP reconstruction of dynamic data. (a) Reconstructed image. (b) Detail of the vessels. (c) Zoomed detail of the background structures.

Figure 5 . 23 :

 523 Figure 5.23: Convergence curves (RMSD values) for the numerical phantom. Reconstructions: CS-piccs (black dashed line), CS-sbs+t (blue line), CS-wm+tv (red line).

Fig. 5 . 4 :

 54 25(a) and the post-processing strategy in Fig.5.25(b) are Parameter settings for CS reconstructions with a prior deformation field.

Figure 5 . 24 :

 524 Figure 5.24: Vector image of the deformation field (a) estimated with DRAMMS algorithm between phase 1 (b) and phase 2 (c). Largest displacements are represented by yellow or orange arrows.

Figure 5 .

 5 [START_REF] Chen | Prior Image Constrained Compressed Sensing (PICCS)[END_REF] shows the deformation field along z-axis that was obtained using DROP within a RoI centered around the catheter. Results on clinical data are analyzed in Fig.5.27 where standard clinical (i.e. static) reconstruction shown in Fig.5.27(a) is compared to CS A+T reconstruction of

Figure 5 . 25 :

 525 Figure 5.25: Reconstruction results with prior deformation field on numerical simulations. (a) Static reconstruction. (b) Motion correction obtained with post-processing deformation field. (c) CS A+T reconstruction.

Figure 5 . 26 :

 526 Figure 5.26: Deformation field along z-axis for the catheter motion study.

Figure 5 . 27 :

 527 Figure 5.27: Axial slice of the catheter motion study. (a) Static reconstruction. (b)-(d) CS A+T reconstruction for different phases.

  , A. Reshef, C. Riddell, Y. Trousset, A. Tenenhaus, E. Lahalle, G. Fleury, N. Paragios. Nonlinear diffusion constraints for reconstructing subsampled rotational angiography data. Proceedings of the 12th International Meeting on Fully Three-dimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully3D'2013), 4 pp., Lake Tahoe, California, 16-21 June 2013. H. Langet, C. Riddell, Y. Trousset, A. Tenenhaus, E. Lahalle, G. Fleury, N. Paragios. Compressed Sensing Based 3D Tomographic Reconstruction for Rotational Angiography. Medical Image Computing and Computer Assisted Intervention (MICCAI'2011), 8 pp., Toronto, Canada, 18-22 September 2011. • H. Langet, C. Riddell, Y. Trousset, E. Lahalle, A. Tenenhaus, G. Fleury, N. Paragios. Sparsity constraints and dedicated acquisition protocols for improved Digital Subtraction Rotational Angiography. Proceedings of the 11th International Meeting on Fully Three-dimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully3D'2011), 4 pp., Potsdam, Germany, 11-15 July 2011. Langet, C. Riddell, A. Tenenhaus, Y. Trousset, E. Lahalle, G. Fleury, N. Paragios. Threedimensional tomographic reconstruction with ℓ 1 -minimization strategy for rotational angiography. 20th International Conference on Computational Statistics (COMPSTAT'2012), Limassol, Cyprus, 27-31 August 2012.

	Invited talks
	• H.
	• H. Langet, C. Riddell, Y. Trousset, A. Tenenhaus, E. Lahalle, G. Fleury, N. Paragios. Compressed Sensing Dynamic Reconstruction in Rotational Angiography. Medical Image Computing
	and Computer Assisted Intervention (MICCAI'2012), 8 pp., Nice, France, 1-5 October 2012.

• H. Langet, C. Riddell, Y. Trousset, A. Tenenhaus, E. Lahalle, G. Fleury, N. Paragios. Compressed sensing subtracted rotational angiography with multiple sparse penalty. Proceedings of the 9th International Symposium on Biomedical Imaging (ISBI'2012), 8 pp., Barcelona, Spain, 2-5 May 2012. • Patents • H. Langet, C. Riddell, for General Electric Company. Tomographic reconstruction of an object in motion, EP2469471A1 / US20120163532, 27 June 2012.

  2.1.

	Material	Linear attenuation coefficients (in cm -1 )	Hounsfield unit
	Titanium	1.23	6200
	Cortical bone Whole blood Lung Muscle Brain (grey and white matter) Water Adipose tissue Air	3.56 •10 -1 1.80 •10 -1 1.78 •10 -1 1.78 •10 -1 1.77 •10 -1 1.71 •10 -1 1.60 •10 -1 1.86 •10 -4	1080 50 40 40 35 0 -65 -1000
	Table 2.1: Linear attenuation coefficients at 100 keV (source: NIST tables of X-ray mass attenuation
	coefficient) and corresponding Hounsfield units for different materials.

  • /s with Innova and IGS systems) for acquisition with injection of a contrast agent. The acquisition at 30 fps and 40 • /s thus delivers a scan of about 150 views with an acquisition time of about 5 seconds since the gantry can only rotate over 200 • . On the other hand, the rotation of the gantry is kept slow for non-angiographic exam (10 • /s with Innova systems, 16

• /s with IGS systems) in order to maximize the amount of views within the scan. Acquisitions at 30 fps and 10 • /s and acquisitions at 50 fps and 16 • /s thus deliver scans of about 600 views. Both Innova and IGS systems are equipped with an intermediate rotation speed (20 • /s with Innova systems, 28

  p 11 p 12 p 13 p 14 p 21 p 22 p 23 p 24 p 31 p 32 p 33 p 34

Table 3 .

 3

	Acquisition time	20 s -10 s -5 s	12.5 s -7 s -5 s
	Amount of views	600 -300 -150	600 -350 -200
	Angular coverage	200 •	200 •

• /s -40 • /s 16 • /s -28 • /s -40 • /s 1: Performance of C-arm systems.

Table 4 .

 4 2: Parameters of the sparsity constraint and RMSD values.

	Experiment Algorithm Regularization	Homotopy	Final RMSD value
	n • 1	iFBP+			0.074
	n • 2	CS-tv	λ tv = 10		0.059
	n • 3	CS-tv	λ tv = 1		0.065
	n • 4	CS-w	λ w = 10		0.098
	n • 5	CS-w	λ w = 1		0.065
	n • 6 n • 7 n • 8 n • 9	CS-sbs CS-tv CS-tv CS-w	λ λ (s) sbs = 3000 to 0 (s) tv = 3000 to 0 λ (s) tv = 3000 to 1 λ (s) w = 3000 to 0	0.066 0.058 0.057 0.069
	n • 10	CS-w	λ	(s) w = 3000 to 1	0.068

• 7 and n • 9) and TV and wavelet penalties (experiments n • 8 and n • 10) reconstructions. The quality of these reconstructions is compared in Fig. 4.10b and Fig. 4.10c for TV, and in Fig. 4.10d and Fig. 4.10e for wavelet. These figures indicate that homotopy reconstructions are superior to regularized

Table 5 .

 5 2: Results in parallel geometry on a numerical phantom with simulated motions. RMSD values.

		No motion	Pulsatile motion	Successive motion
		Static	Static	CS-sbs+t Static CS-sbs+t
	d 1	35 HU	187 HU	68 HU	91 HU	85 HU
	d 2	44 HU	159 HU	71 HU	83 HU	83 HU
	d 3	49 HU	144 HU	73 HU	81 HU	80 HU

Hence, the term 'tomography' which is formed from the Greek roots 'tomos' (slice) and 'graphein' (write).

Anti-scatter collimator is also referred to as 'grid' or 'septa'.

Arteriovenous malformations refer to abnormal connections between arteries and veins

Coil herniation is a complication in which the coil goes out of the aneurysm into the parent artery.

Noticing that the coefficient of a given scale is likely to be greater than the coefficient of the same orientation in the same spatial location at finer scales, if a coefficient is set to zero, its descendens at finer scales are also set to zero.

RMSD definition is given in Sec.3.1.4.3.

RMSD definition is given in Sec.

3.1.4.3.

Since DRAMMS algorithm is a 3D registration algorithm, the 2D numerical phantom was replicated 20 times to form a thin volume.

Remerciements

Abstract:

Medical imaging has known great advances over the past decades to become a powerful tool for the clinical practice. It has led to the tremendous growth of interventional radiology, in which medical devices are inserted and manipulated under image guidance through the vascular system to the pathology location and then used to deliver the therapy. In these minimally-invasive procedures, X-ray guidance is carried out with C-arm systems through two-dimensional real-time projective low-dose images. More recently, three-dimensional visualization via tomographic acquisition has also become available.

This work tackles tomographic reconstruction in the aforementioned context. More specifically, it deals with the correction of motion artifacts that originate from the temporal variations of the contrast-enhanced vessels and thus tackles a central aspect of tomography: data (angular) sampling. The compressed sensing theory identifies conditions under which subsampled data can be recovered through the minimization of a least-square data fidelity term combined with sparse constraints. Relying on this theory, an original reconstruction framework is proposed based on iterative filtered backprojection, proximal splitting, ℓ 1 -minimization and homotopy. This framework is derived for integrating several spatial and temporal penalties. Such a strategy is shown to outperform the analytical filtered backprojection algorithm that is used in the current clinical practice by reducing motion and sampling artifacts in well-identified clinical cases, with focus on cerebral and abdominal imaging. The obtained results emphasize one of the key contributions of this work, that is the importance of homotopy in addition to regularization, to provide much needed image quality improvement in the suggested domain of applicability.

Keywords: Tomographic reconstruction, iterative methods, intraarterial angiography, C-arm systems, compressed sensing, ℓ 1 -penalization, proximal splitting, homotopy, registration

Échantillonnage et reconstruction de mouvement en radiologie interventionnelle tridimensionnelle

Résumé :

La pratique clinique a été profondément transformée par l'explosion technologique, ces dernières décades, des techniques d'imagerie médicale. L'expansion de la radiologie interventionnelle a ainsi rendu possible des procédures dites "minimalement invasives" au cours desquelles la thérapie est délivrée directement au niveau de la région pathologique via des micro-outils guidés par imagerie à travers le système vasculaire. Des systèmes dits "C-arm", générant une imagerie rayons X planaire temps-réelle en faible dose, sont utilisés pour le guidage. Ils ont offert plus récemment la possibilité d'une visualisation tridimensionnelle par le biais d'acquisitions tomographiques.

C'est dans ce contexte de reconstruction tomographique que s'inscrivent ces travaux de thèse. Ils s'attèlent en particulier à corriger les artefacts de mouvement dus aux variations temporelles des vaisseaux injectés et se concentrent sur un aspect central de la tomographie, à savoir l'échantillonnage angulaire. La théorie du compressed sensing identifie les conditions sous lesquelles des données sous-échantillonnées peuvent tre reconstruites en minimisant une fonctionnelle qui combine un terme de fidélité quadratique et une contrainte parcimonieuse. S'appuyant sur cette théorie, un formalisme original de reconstruction est proposé : il repose sur la rétroprojection filtrée itérative, les algorithmes proximaux, la minimisation de normes ℓ 1 et l'homotopie.

Ce formalisme est ensuite dérivé pour intégrer différentes contraintes spatiales et temporelles. Une telle stratégie s'avère plus performante que la rétroprojection filtrée analytique utilisée dans la pratique clinique, permettant la réduction d'artefacts de mouvement et d'échantillonnage dans des cas cliniques bien identifiés de l'imagerie cérébrale et abdominale. Les résultats obtenus soulignent l'une des principales contributions de ce travail, à savoir : l'importance de l'homotopie, en supplément de la régularisation, pour améliorer la qualité image, un gain indispensable dans le domaine d'applicabilité suggéré.

Mots clés : Reconstruction tomographique, méthodes itératives, angiographie intra-artérielle, systèmes C-arm, compressed sensing, pénalisation ℓ 1 , algorithmes proximaux, homotopie, recalage