A B S T R A C T

En imagerie par transmission de rayons X, les mesures effectuées relèvent d'un opérateur intégral : la transformée de Radon en géométrie parallèle et la transformée conique en géométrie divergente. La caractérisation de l'image de ces opérateurs par des conditions de rang permet de quantifier l'incohérence des données mesurées par rapport au modèle intégral. Le premier volet de ce travail de thèse étudie les conditions de rang en géométrie conique : nous proposons de nouvelles conditions pour une trajectoire planaire et mettons à jour des relations non triviales entre des conditions 2D et des conditions basées sur le théorème de Grangeat. Le second volet porte sur l'auto-étalonnage géométrique des systèmes tomographiques à géométrie conique. L'analyse des conditions de rang couplée au modèle géométrique des projections radiographiques permet la détermination de la géométrie d'acquisition du système.

In X-ray transmission imaging, the collected measurements correspond to an integral operator: the Radon transform in parallel geometry and the divergent beam transform in divergent geometry. The range of these operators is characterized by conditions, which help to quantify the consistency of the measured data with the forward integral model. The first pillar of this PhD work studies range conditions in cone-beam acquisition geometry: we derive new conditions for a planar trajectory and establish a new relation between 2D fanbeam conditions and Grangeat-based conditions. O,e e x , e e y , e e z ). Applies to a pair of projections.

DCS Detector Coordinate System. (O d , u, v ). The coordinate system in the same unit as the PCS, with origin at the centre of the detector and axes in the same direction as the x-and y-axes of the OCS.

3DCS 3D Detector Coordinate System. (O d , u, v , w ). It is a DCS turned into a 3D Coordinate system by adding a unit vector orthogonal to the plane of the detector.

RHS Right hand-side (of an equation).

LHS Left hand-side (of an equation).

xiii I N T R O D U C T I O N

Radiology is the medical science that uses the three main modalities of medical imaging: radiation, magnetic resonance and ultra-sound. This work's main interest is on radiation imaging. The radiation imaging field itself is divided into two main categories: transmission imaging and emission imaging. In emission imaging, a radioactive product is injected into the patient. The interaction of the product with the potentially abnormal cells of the patient's organ leads to the emission of photons. Here, we focus on transmission imaging (though some theoretical concepts presented in this work remain valid in emission imaging). In such devices, a beam of X-rays is emitted from a X-ray generator and traverses the body of a patient. As photons go through the tissues of the body, they interact with matter and are attenuated. The outgoing attenuated X-rays are collected on a detector. One single image is called a radiograph. When the process is repeated several times at different view angles, the attenuation map of the patient's organ can be numerically reconstructed from the set of all radiographs. This reconstruction process is called computerized tomography (CT). The basics of CT will be presented in Part I, Chapter I-1.

There are several types of CT devices. The most widely used is the usual CT scanner (see Figure 1a). A CT scan is installed in a dedicated room. Such scanners are exclusively used for diagnostic. They are mechanically very robust. The gantry rotates at a fairly high frequency (up to a few rotations per second) with very high stability. Cormack and Hounsfield, who built the first CT scanner in 1971 1 were granted the Nobel Prize in medicine for this 1 achievement. In parallel, radiological systems have been designed to enter the Operating Room (OR) and help the clinician not only for the diagnostic but during the treatment. This field is called interventional radiology. The first devices in this direction were mobile C-arms (illustrated in Figure 1b), which were only capable of conventional radiography or fluoroscopy (continuous illumination and visualization of the patient). Thanks to this new technology, diseases which required an open surgery could be treated in a mini invasive way, by the use of image guided catheters. With progress in robotics, manufacturers developed robotic C-arms (see Figure 1c) with 3D CT capability and hospitals have implemented dedicated IR operating rooms. The 3D imaging capabilities of robotic C-arms were made possible by many factors (large flat panel detectors, new reconstruction algorithms, etc.) but one of the challenges was to be able to accurately know the acquisition geometry of each shot. This problem of estimating the geometry is called calibration. It is explained in Part I, Chapter I-2 and a review of existing methods is presented in Part I, Chapter I-4.

The reason why mobile C-arms were not equipped with 3D CT feature is mostly related to the calibration problem. Robotic C-arms or CT scanners are mechanically designed so that the deformations of the system are negligible. A geometric calibration is only necessary periodically (about every six months) to account for the mechanical wear of parts. This calibration relies on the crucial assumption that the acquisition geometry is reproducible over time. On the contrary, mobile C-arms are subject to unpredictable and non-reproducible vibrations and deformations due to their light-weight design. Hence the need for calibration procedures which can be carried at each production scan. As explained in Chapter I-4, one possible way of achieving this is to place fiducial markers in the scene with known respective positions. To avoid cluttering the medical scene, the ideal solution would rely only on the analysis of the projection data. Data Consistency Conditions (DCC) are mathematical conditions which quantify the gap between the measurements (the projection data) and the theoretical forward line integral model. A poor estimation of the geometry leads to a poor consistency of the data. This observation paves the way to on-line calibration. This PhD work is twofold. On one hand, we deeply studied the mathematical aspects of DCC, derived new cone-beam DCC and established the unknown relation between two known sets of DCC. On the other hand, we applied DCC to the calibration of a micro-CT system. This calibration procedure is -hopefully -a first step towards pure on-line calibration of a mobile C-arm.

The manuscript is organized as follows. Part I presents all the scientific background which the reader should know before reading the contributions. Chapter I-1 presents the mathematics of tomography. In Chapter I-2, we describe the acquisition geometry, introduce the calibration problem and define the projection matrix. Chapters I-4 and I-3 present a review of the literature on DCC and geometric calibration respectively. In Part II, we gather the three contributions of this work. Chapter II-1 presents new conebeam DCC. This contribution partially relies on the theory of distributions whose basics are presented in Appendix 1.1. Chapter II-2 presents a new on-line DCC-based calibration method for a circular cone-beam turntable CT system. Finally, Chapter II-3 establishes the relation between the fanbeam DCC used in our second contribution and Grangeat-based DCC.

Part I S C I E N T I F I C B A C K G R O U N D

This part gathers all the scientific knowledge that a nonspecialist reader needs for a smooth reading of the contributions. It also places the contributions in their scientific context and defines the notation.

I1 T O M O G R A P H Y : A N O V E R V I E W

The aim of tomography is to reconstruct a function of interest from data (called projections), whose formation process is known and precisely modelled. This model is called the forward model. Such a problem is called an inverse problem. In this chapter, we introduce the various forward models that are used in this work. We also introduce their dual operators and briefly describe some basic reconstruction methods. The last section gives a proof of the Grangeat theorem. This theorem plays a central role in the work and its proof gently introduces the mathematical methods involved in the contributions.

T H E B E E R -L A M B E R T L A W

In transmission imaging, X-rays are emitted by an external point X-ray source, traverse the patient modelled by its density map µ, and are collected on a 2D flat-panel detector. The image output is a radiograph. The aim of tomography is to recover the density map µ from radiographs. The density map µ may be defined on R 2 (single slice reconstruction) or R 3 . The correct assumption to make on µ would be: bounded with compact support. Many of the mathematical results presented in this work would remain valid with such µ. Nevertheless, the presentation would be slightly more tedious. So we simply assume that µ ∈ D(R), the space of compactly supported smooth functions 2 , which itself is contained in S (R n ) 2 Compactly supported bounded functions are in L 1

(R n ) and D(R n ) is dense in L 1 (R n )
(n ∈ {2, 3}). When an object µ is irradiated with X-rays with initial intensity I 0 , the outgoing X-ray (which is in practice being collected by a detector) has intensity I . The Beer-Lambert law links the measured intensity I to the density map µ as

I = I 0 exp - L µ( x) d x , (1)
where photons have travelled along the line L (from the emitting source to the detector's pixel). Since I 0 is a characteristic feature of the scanner (X-ray source and detector) and does not depend on the object µ, the reconstruction problem boils down to recovering µ from integral of µ along lines. Given a line L, the collected data will be modelled in this work by the negative log transform of the physically collected data:

-log I /I 0 = L µ ( x) d x. ( 2
)
Different scanning geometries correspond to different ways of collecting line integrals. This is discussed in the next section. Unfortunately, the line integral model is ideal in the sense that the collected data on real devices are in practice always degraded with systematic effects such as beam hardening or scattering amongst others, which introduce artefacts in the reconstructed image. This work only addresses one of these problems: the geometric calibration. The general idea in this work is to carry out these two steps:

1. quantify the gap between the collected data and the ideal model (see Data Consistency Conditions in Chapter I-3), 2. minimize this gap by adjusting the geometric parameters which are responsible for this deviation from the ideal model (see Calibration in Chapter I-4).

S O M E C O O R D I N AT E C O N V E N T I O N S

Throughout this document, we will often switch between different coordinate systems. All coordinate systems are illustrated in Figure 2. See also [START_REF] Siddon | Solution to treatment planning problems using coordinate transformations[END_REF].

• The Patient Coordinate System (PCS) is the generic 3D coordinate system. Its origin is O and the unit vectors are ( e x , e y , e z ). The associated coordinates are denoted (x, y, z).

• The Detector Coordinate System (DCS). Given a detector (e. g. a flatpanel pixel grid), the origin of the DCS is O d , the centre of the rectangular panel. The unit vectors are ( e u , e v ). They coincide with the detector rows and columns respectively. The associated coordinates are usually denoted (u, v).

• The 3D Detector Coordinate System (3DCS) is an orthonormal extension of 2D DCS to a 3D right-handed coordinate system. More precisely, if (O d , e u , e v ) is a DCS. We set e w = e u × e v ). The 3DCS is (O d , e u , e v , e w ).

• The Oriented Coordinate System (OCS). Its origin is O and the unit vectors are ( e o x , e o y , e o z ). The associated coordinates are (x o , y o , z o ). In general terms, it is a rotated version of the PCS. It is associated to a rotation matrix R such that (x o , y o , z o ) T = R(x, y, z) T . In practice, it is related to a DCS in the following way: a 3DCS and an OCS have the exact same triple of unit vectors. They only differ by the origin.

F O R W A R D O P E R AT O R S

We introduce several integral transforms which model the variety of possible scanning geometries. The Radon transform maps a density function µ to the set of its integrals over hyperplanes.

Definition 1 (The Radon transform). Let µ ∈ S (R n ). The Radon transform of µ is defined by

Rµ( θ, s) = x• θ=s µ( x) d x = R n δ x • θ -s µ ( x) d x, θ ∈ S n-1 , s ∈ R. (3)
In the last integral, δ denotes the Dirac δ distribution (see Appendix 1.1 for an introduction to δ and some of its properties). There is alternative notation for the integral in Equation 3,

Rµ( θ, s) = θ ⊥ µ s θ + y d y, θ ∈ S n-1 , s ∈ R, (4)
and another, in the particular 2D case:

R 2 µ( θ, s) = R µ s θ + l ζ ⊥ dl , θ, ζ ∈ S 1 such that , θ • ζ = 0. ( 5
)
Note that the variables of the Radon transform identify one hyperplane:

s θ φ Figure 4
θ is the normal unit vector to the hyperplane and s is its signed distance to the origin (see Figure 4). Note that this identification is not unique.

We denote H θ,s such hyperplanes. In 2D, a hyperplane is a line and in 3D, it is a plane. Note also that in the 2D case, the unit vector θ can be easily parametrized by the angle that the ray makes with the x-axis. Let φ ∈ [0, 2π[ denote that angle. Then, we have θ = (-sin φ, cos φ) and we will use interchangeably θ and φ as the first parameter. Finally, in the 2D case, the Radon transform describes the physical reality of X-ray measurements. In 3D, the Radon transform does not model the X-ray acquisition. We need to introduce the X-ray transform.

Definition 2 (The X-ray transform). Let µ ∈ S (R 3 ). The X-ray transform of µ is defined by

X µ( θ, y) = R µ y + l θ dl , θ ∈ S n-1 , y ∈ θ ⊥ , ( 6
)
where θ ⊥ = y ∈ R n / y • θ = 0 .

These two transforms -R and X -are relevant with parallel acquisition geometries: within one projection, all rays are parallel. In the case of X-rays diverging from a point source, we define the divergent transform: Definition 3 (The divergent-beam transform). Let µ ∈ S (R n ). Let ( s λ ) λ∈Λ⊂R denote the 1D trajectory of a X-ray source. The divergent-beam transform is defined by:

Dµ λ, α = ∞ 0 µ s λ + t α dt , λ ∈ Λ, α ∈ S n-1 . ( 7
)
We suppose that for all λ ∈ Λ, the source is outside the convex hull of the support of µ, so that the integration could equivalently be over the entire real line R. For the applications, the circular trajectory is of particular importance. Without loss of generality, we consider the circular trajectory centred on the origin of the PCS and in the plane y = 0 in the 3D case. The parameter λ thus reduces to the polar angle of the source with respect to the x-axis in the trajectory plane. In the 3D case, referred to as cone-beam, the unit vector α is parametrized using spherical coordinates. Such parametrization will be seen extensively in the next sections. In the 2D case, referred to as fanbeam, the unit vector α is parametrized with the angle made by the ray with respect to a reference axis. Two common situations are:

• The reference axis is the x-axis of the PCS as for the 2D Radon transform. The angle is denoted φ. .

• The reference axis is the central ray (joining the source position and the origin of the PCS). The angle is denoted γ. See Figure 5.

S O M E N O TAT I O N T I P S

We consistently use the first variable to identify one single projection. The second variable identifies one ray or hyperplane in that projection. When considering one single projection, we may use some specific notation. For example, p is usually used to denote the Radon transform of a function µ. The projection at fixed θ is denoted p θ . For divergent data, we use the letter g . One projection along the trajectory is denoted g λ . We have

p θ (s) = p θ, s = Rµ θ, s (8)
g λ α = g λ, α = Dµ λ, α . (9)
We call re-binning the operation that maps divergent-beam data to parallel-beam data. More specifically, in the two particular ray parametrizations in fanbeam geometry described above, the re-binning operations are given by g λ, φ = p φ, R sin λφ [START_REF] Clackdoyle | Tomographic Reconstruction in the 21st Century[END_REF] g λ, γ = p λγ, R sin γ , [START_REF] Clackdoyle | Fanbeam data Consistency Conditions for applications to motion detection[END_REF] where g and p are divergent-beam and parallel-beam data respectively.

1.4 T H E F O U R I E R S L I C E T H E O R E M
The reader is referred to the textbook of Natterer [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF] for a complete presentation of these operators and their properties. We just mention the obvious fact that the transforms are linear operators (because the integral is linear). They are defined in terms of direction of the integration line, regardless of any physical detector. In the applications, a detector will be considered and the second variable will be substituted with the detector's variables.

1.4 T H E F O U R I E R S L I C E T H E O R E M
In view of this major (though very simple) theorem, we will need the Fourier transform, which is extensively used in this work. For our purpose, it will be sufficient to introduce the Fourier transform on the Schwartz spaces S (R n ) of C ∞ rapidly decreasing functions 3 , which conveniently

∈ Z n + , l ∈ Z n + , x k ∂ |l | f ∂x l → 0 as |x| → ∞
maps onto itself. See e. g. [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I[END_REF] for a comprehensive treatment of the Fourier transform

. Let f ∈ S (R n ). Since f is integrable on R n , the functions f ξ = R n f ( x) e -2i π ξ• x d x (12) f ( x) = R n f ξ e 2i π ξ• x d ξ (13)
are well defined. A key theorem in Fourier analysis (which we do not prove) is

Theorem 1 (Fourier Inversion formula). For f ∈ S (R n ), f ξ = f ξ ∀ ξ ∈ R n (14) f ( x) = f ( x) ∀ ξ ∈ R n . ( 15
) ( 16
)
We can now define the Fourier transform and its inverse Definition 4 (The Fourier transform). Let f ∈ S (R n ). The Fourier transform of f is defined over R n by

F f ξ = f ξ = R n f ( x) e -2i π ξ• x d x. ( 17
)
Its inverse is

F -1 f ( x) = f ( x) = R n f ξ e 2i π ξ• x d ξ. ( 18
)
We gather in the next Proposition, which we do not prove, the properties of the Fourier transform which we will need in the sequel. For a multi-index

k ∈ Z n + = (k 1 , ..., k n ), |k| = k 1 +...+k n , x k = x k 1 1 ...x k n n and f (k) = ∂ k 1 1 ...∂ k n n f . The standard convolution product is denoted f * g and defined as f * g y = R n f ( y -x)g ( x) d x = R n f ( x)g ( y -x) d x. Proposition 2 (Properties of the Fourier transform). Let f , g ∈ S (R n ). F f (k) (ξ) = 2i π) |k| ξ k F f ξ (19) F x k f ξ = (2i π) |k| F f (k) ξ (20) f * g = f ĝ (21) f g = f * ĝ . ( 22
)
Most of the time in this work, when applied to projection functions, the Fourier transform acts on the second variable. The next theorem is a cornerstone in image reconstruction.

Theorem 3 (The Fourier slice theorem). Let θ ∈ S n-1 and let p θ (s) denote the Radon transform of a function µ in the direction θ. Then:

p θ (σ) = μ σ θ (23)
Proof. By definition of the Fourier transform

p θ (σ) = R p θ (s) e -2i πsσ ds (24) = R θ ⊥ µ s θ + y e -2i πsσ d y ds (25) = R n µ ( x) e -2i π x• θ σ d x (26) = R n µ ( x) e -2i π x• σ θ d x, (27)
where we applied the unit Jacobian change of variables x = s θ + y and noticed that s = x • θ.

The Fourier slice theorem is the basis of many reconstruction techniques. In particular, the Filtered Backprojection (FBP) is directly derived from this theorem (see Section 1.6). Before moving onto reconstruction, it is necessary to introduce the backprojection operators.

T H E D

U A L O P E R AT O R S : T H E B A C K P R O J E C T I O N
We quickly give (with no proof, see [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF]) the dual operators of the forward operators defined in the previous section.

Proposition 4. For µ ∈ S (R n ), the dual operators of R, X and D, denoted R * , X * and D * respectively are for g ∈ S (S i -1 × R), R * g ( x) = S i -1 g θ, x • θ d θ ( 28
)
for g ∈ S (S 2 × R 2 ), X * g ( x) = S 2 g θ, x -x • θ θ d θ ( 29
)
for g ∈ S (Λ × S 2 ), D * g ( x) = Λ g λ, x -s λ x -s λ dλ x -s λ n-1 . (30)
These operators are called backprojection operators because they correspond to the smearing of the projection values back along the lines (or planes) of integration. After the backprojection is complete, one voxel value will be the sum of all projection values stemming from lines passing through that voxel. The backprojection is not an inverse operator. It will be seen in the next section that projections need to be pre-processed before being back projected. For now, we simply illustrate the fact that the pure backprojection of the projection data leads to a blurred reconstructed image, as shown in Figure 6.

T H E B A S I C S O F R E C O N S T R U C T I O N

In this section, we gather the basic facts about reconstruction. Especially, we derive the so-called filtered backprojection formulas in 2D and 3D, directly from the Fourier Slice theorem. Again, the reader is referred to [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF] for a comprehensive treatment of the inversion of the Radon transform and a family of inversion formulas based on the Riesz potential. For computer implementation of these reconstruction methods, a useful reference would be the book of Kak and Slaney [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. The textbook of Deans [START_REF] Deans | The Radon Transform and some of its applications[END_REF] is also an excellent introduction to the topic.

Inversion of the Radon transform

We begin with

Theorem 5 (Filtered backprojection). Let p ∈ S (S n-1 × R). If p = Rµ then µ( x) = 1 2 S n-1 +∞ -∞ |σ| n-1 p θ (σ) e 2i π x•σ θ dσ d θ. ( 31
)
Proof. Starting from the inverse Fourier transform, changing from Cartesian to polar coordinates and applying the Fourier Slice Theorem, we have

µ( x) = R n μ ξ e 2i π x• ξ dξ (32) = S n-1 +∞ 0 |σ| n-1 μ σ θ e 2i π x•σ θ dσ d θ (33) = S n-1 +∞ 0 |σ| n-1 p θ (σ) e 2i π x•σ θ dσ d θ, ( 34
)
where p θ denote the (1D) Fourier transform of the projection at fixed direction θ with respect to the spatial variable s. Note that the use of the modulus in the inner integral is not necessary since σ > 0 but is intended in view of computing the integral:

S n-1 0 -∞ |σ| n-1 p θ (σ) e 2i π x•σ θ dσ d θ ( 35
)
By substituting (σ, θ) with (-σ, -θ) in this integral, and using the evenness of the Radon transform (p -θ (-σ) = p θ (σ)) 4 , we get that this integral is equal to the integral in the RHS of Equation 34. By adding the two terms, we finally obtain the desired result.

Equation 31 gives a complete method for the recovery of the object function µ in two steps:

1. Modify the projections by multiplying their Fourier transforms by |σ| n-1 , 2. Back-project the modified projections.

It is necessary to further study the modification operations in the 2D and the 3D cases. In the 3D case (n = 3), Equation 31 can be re-written

µ( x) = - 1 8π 2 S 2 +∞ -∞ (2i πσ) 2 p θ (σ) e 2i π x•σ θ dσ d θ (36) = - 1 8π 2 S 2 ∂ 2 ∂s 2 p θ (s) s= x• θ d θ, ( 37
)
so that the modification operation is actually a differentiation with respect to the second variable.

To examine the 2D case, we need the following Definition 5 (The Hilbert transform). For f ∈ S (R), the Hilbert transform of f , denoted H f , is defined by

H f (s) = +∞ -∞ f (s ) π(s -s ) ds (38)
Let us denote h(s) = 1/(πs) the convolution kernel in the above definition. It can be shown that the Hilbert transform in the spatial domain is equivalent to multiplying in the Fourier domain by (-i sign(σ)), where the sign function is defined by

sign(σ) =      σ/|σ| if σ = 0 0 otherwise. (39)
We are now ready to address the 2D FBP. Starting again from Equation 31 with n = 2, we have

µ( x) = 1 4π S 1 +∞ -∞ (2i πσ) -i sign(σ) p θ (σ) e 2i πσ x• θ dσ d θ, ( 40
) = 1 4π S 1 ∂ ∂s H p θ x • θ d θ. ( 41
)
We can now pinpoint the critical difference between the 2D and the 3D cases. In the 3D case, the filtering step is a simple differentiation, which is a local operation (it only involves a small neighbourhood around the value of interest). Consequently, reconstructing the function at a point x only requires the knowledge of the projection data over planes traversing a small neighbourhood of x. On the other hand, in the 2D case, the Hilbert transform is a convolution of the projection data with a convolution kernel h which is nowhere zero. Hence, the computation of H p θ requires all the projection data. Truncated data cannot be handled using a direct implementation of FBP. The problem of reconstructing an image from truncated data is of major interest. It has long been considered that incomplete data necessarily introduced artefacts in the whole reconstructed image. A breakthrough was proposed in [START_REF] Noo | Image reconstruction from fan-beam projections on less than a short scan[END_REF]: a region of interest could be exactly reconstructed from incomplete data. See [START_REF] Clackdoyle | Tomographic Reconstruction in the 21st Century[END_REF] for a historical review. Note that we have described the modification of the projection in the 2D case as a convolution with the Hilbert kernel followed by a differentiation. We could have condensed both into one filtering (convolution) operation by using the so-called ramp filter, defined in the Fourier space by ρ(s) = |σ| e 2i πσs dσ.

(

) 42
The modification step

+∞ -∞ |σ| p θ (σ) e 2i π x•σ θ dσ ( 43
)
can then be rewritten in the form of a convolution ρ * p θ x • θ , hence the name of filtered backprojection.

Fan-beam filtered backprojection and the FDK algorithm

In this section, we focus on the two-dimensional case and derive an FBP formula for fanbeam projections acquired along a circular trajectory. The general idea is to apply the FBP inversion formula of Equation 31 -valid for parallel projections -to the fanbeam data after a change from fanbeam variables to the Radon space variables. We assume a circular trajectory of the source, with radius r and let λ denote the angle made by the source with respect to the y-axis. Therefore, s λ = (-sin λ, cos λ). For convenience of the derivation, the 1D detector is aligned 6 of the PCS. A pixel's position on the detector is identified with its abscissa u. If α(u) denotes the unit vector originating at s λ in the direction of the detector position u (see Figure 7), we have

g (λ, u) = +∞ 0 µ( s λ + t α(u)) dt . ( 44
)
In the 2D Radon space, we identify the direction of a projection θ with the angle φ that θ makes with the y-axis. We have the following re-binning formula from fanbeam coordinates to Radon space coordinates:

s = r u r 2 + u 2 (45) θ = λ + tan -1 u r ( 46
)
Applying this change of variables in Equation 31 yields the following inversion formula:

Theorem 6 (Fan-beam FBP). µ( x) = 2π 0 1 M 2 +∞ -∞ g (λ, α(u)) r r 2 + u 2 ρ (u x -u) du dλ, ( 47
)
where M is the magnification factor M = r -x• w r , u x = x• u M and ρ the ramp filter defined in Equation 42.

The fully detailed derivation of the previous fanbeam FBP formula can be found in [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. We simply outline here the physical meaning of each term in the formula. First, the projection data are weighted with the cosine of the incidence angle of the ray (the r / r 2 + u 2 term). The weighted projections are then ramp-filtered. Finally the filtered projections are backprojected with a 1/M 2 weight.

From fanbeam FBP to FDK

Though this thesis's main interest is on 3D, the 2D fanbeam FBP reconstruction formula is of major importance because it is the basis for the most popular analytic 7 reconstruction algorithm for 3D reconstruction from cone- 7 The qualifier analytic here means that it is an explicit inversion formula. It does not mean that the algorithm is exact.

beam projections: the Feldkamp-Davis-Kress (FDK) algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF].

A concise presentation of the FDK algorithm can be found in [START_REF] Rodet | The cone-beam algorithm of Feldkamp, Davis, and Kress preserves oblique line integrals[END_REF], from which we extract the following outline. The cone-beam data are acquired along a circular trajectory in the y = 0 plane (the axis of rotation is (O y)), with radius r . At each projection λ, the 2D flat-panel detector is aligned and placed at the origin O. Heuristically, the FDK algorithm considers detector rows which are not in the plane of the trajectory as 1D fanbeam projections in the plane of a another (virtual) circular trajectory. Of course, this has to be compensated for in the treatment of the corresponding line integral values. Using the notation of PCS and 3DCS of the detector, the algorithm numerically implements the following formula: Definition 6 (The FDK reconstruction formula).

µ F DK ( x) = 1 2 2π 0 1 M 2 +∞ -∞ g λ (u, v) r r 2 + u 2 + v2 ρ( ū -u) du dλ, ( 48
)
where ( ū, v) are the coordinates of the intersection of the detector with the ray from s λ passing through x. Explicitly, we have:

ū = r x • u r -x • w (49) v = r x • v r -x • w (50)
The FDK algorithm is very similar to the fanbeam FBP algorithm. We emphasize the main steps:

• Weight the projections with the cosine of the incidence angle,

• Ramp-filter each detector row (see Theorem 6 for the definition of the ramp filter),

• Weight the filtered projections with the square of the inverse of the magnification factor,

• Back-project the obtained filtered and weighted projection data.

We emphasize that the FDK algorithm is exact in the plane of the trajectory (this means that if g = Dµ, then µ FDK = µ for y = 0) but only approximates the sought function µ outside this plane. More can be said on the accuracy of FDK, see for example [START_REF] Rodet | The cone-beam algorithm of Feldkamp, Davis, and Kress preserves oblique line integrals[END_REF].

Cone-beam reconstruction

The FDK algorithm described in the previous section falls in the category of cone-beam reconstruction but is essentially a heuristic extension of the 2D FBP algorithm to the 3D case. In this section, we address the challenging question of recovering a 3D object function µ from 2D cone-beam projections. We assume that each 2D projection is not truncated (the support of the object function µ is fully contained in the field-of-view of the scanner) and that the source trajectory is complete in the sense of Tuy [START_REF] Tuy | An Inversion Formula for Cone-Beam Reconstruction[END_REF]: each plane crossing the support of the object function µ intersects nontangentially the trajectory of the source. This condition roughly guarantees that for each point x to be reconstructed, all planes passing through a neighbourhood of x are measured so that the 3D Radon inversion formula of Equation 37 can be applied. A two-orthogonal circles trajectory satisfies Tuy's condition. On the other hand, the single circle trajectory does not, since planes parallel to the trajectory have no intersection with the trajectory. The completeness conditions are not discussed here. We simply mention that most clinical CT scanners (with helical trajectory) satisfy Tuy's condition but C-arm (circular trajectory) do not. Instead of sequentially presenting the various cone-beam inversion methods of Tuy [START_REF] Tuy | An Inversion Formula for Cone-Beam Reconstruction[END_REF], Smith [START_REF] Smith | Image Reconstruction from Cone-Beam Projections: Necessary and Sufficient Conditions and Reconstruction Methods[END_REF] and Grangeat [START_REF] Grangeat | Analyse d'un systeme d'imagerie 3d par reconstruction a partir de radiographies x en geometrie conique[END_REF][START_REF] Grangeat | Mathematical framework of cone beam 3D reconstruction via the first derivative of the radon transform[END_REF], we adopt the unifying framework of [START_REF] Clackdoyle | Overview of reconstruction algorithms for exact cone-beam tomography[END_REF] and [START_REF] Clack | Cone-beam reconstruction by the use of Radon transform intermediate functions[END_REF]. It will give us the opportunity to introduce some intermediate functions, which will play a major role in the contributions of this PhD work. The theoretical framework encompasses the three aforementioned cone-beam inversion methods. All these methods are based on the 3D Radon inversion formula of Equation 37. They differ in the pre-processing of the projection data before backprojection. As usual, the projection data are denoted g (λ, α) where λ ∈ Λ parametrizes the trajectory of the source and α ∈ S 2 is the direction of one ray. The method is two-step:

• For each projection, evaluate an intermediate function G(λ, β) which links the cone-beam data to the Radon transform of the object function µ,

• Back-project a convolved and re-sampled version of the intermediate function.

In view of mathematically deriving the procedure, it is necessary to introduce four real functions (actually, distributions) 8 : the Dirac δ distribution and its derivative (as a distribution) δ , the ramp filter ρ and the Hilbert filter h. Though these functions were already introduced elsewhere, we give a uniform (though not rigorous, see Appendix 1.1) definition in terms of their Fourier transform. For all s ∈ R,

δ(s) = +∞ -∞ e 2i πσs dσ, (51)
δ (s) = 2i π +∞ -∞ σ e 2i πσs dσ, (52)
ρ(s) = +∞ -∞ |σ| e 2i πσs dσ, ( 53
)
h(s) = -i +∞ -∞ sign(σ) e 2i πσs dσ. ( 54
)
These functions have interesting properties: δ and h are homogeneous 9

9 A function f is said homogeneous of degree k ∈ Z if for all a > 0, f (ax) = a k f (x) for all x.
of degree -1 while δ and ρ are homogeneous of degree -2. Also, ρ and δ are even ( f (-x) = f (x), ∀x) while h and δ are odd functions ( f (-x) = f (x), ∀x). The Table 1 summarizes these properties. It is easy to prove (see Appendix ??)that the set of 1D functions, homogeneous of degree -2, form a vector space, that this space has dimension 2, and it is spanned by δ and ρ. We now introduce the intermediate function

G τ (λ, β) = S 2 ε τ α • β g λ, α d α, ( 55
)
where ε τ is a convex combination (0 ≤ τ ≤ 1) of δ and ρ:

ε τ (s) = (1 -τ)δ (s) + τρ(s). ( 56
)
G τ relates the cone-beam projection to the 3D Radon transform of the object function µ through the following derivation, which essentially involves the polar-to-Cartesian change of variables x = s λ +t α, with jacobian

t 2 = x -s λ 2 . G τ (λ, β) = S 2 ε τ α • β +∞ 0 µ s λ + t α dt d α, ( 57
) = R 3 ε τ x -s λ x -s λ • β µ ( x) d x x -s λ 2 , ( 58
) = R 3 ε τ x • β -s λ • β µ ( x) d x, ( 59
) = R 3 +∞ -∞ ε τ s -s λ • β δ x • β -s ds µ ( x) d x, (60)
= +∞ -∞ ε τ s -s λ • β R 3 δ x • β -s µ ( x) d x ds, (61)
= +∞ -∞ (1 -τ)δ + τρ s -s λ • β Rµ β, s ds, ( 62
) = τ ρ * Rµ s λ • β -(1 -τ) δ * Rµ s λ • β (63)
Thanks to the symmetry properties of δ and ρ, the intermediate function takes clearly the form of a filtering of the Radon transform of the object over the plane H β, β• s λ .

Following [START_REF] Clackdoyle | Overview of reconstruction algorithms for exact cone-beam tomography[END_REF], we introduce the function

F τ ( β, s) = G τ (λ, β)
where λ is such that s λ • β = s. The latter condition means that the plane H β,s contains the source position s λ . The reconstruction formula is finally given by

f r ( x) = 1 2(2τ 2 -2τ + 1) S 2 +∞ -∞ ε τ x • β -s F τ β, s ds d β. ( 64
)
Note that the outer integration is not a backprojection since it does not integrate over the source trajectory. The integrand function is a filtered version of the function F τ . Each value of τ provides a reconstruction method. We emphasize here two specific values. For τ = 0, ε 0 = δ and the corresponding reconstruction method is that proposed by Grangeat in [START_REF] Grangeat | Analyse d'un systeme d'imagerie 3d par reconstruction a partir de radiographies x en geometrie conique[END_REF][START_REF] Grangeat | Mathematical framework of cone beam 3D reconstruction via the first derivative of the radon transform[END_REF]. The related intermediate function G 0 will be denoted G g in the sequel. For τ = 1, ε 1 = ρ and the reconstruction method is that of Smith [START_REF] Smith | Image Reconstruction from Cone-Beam Projections: Necessary and Sufficient Conditions and Reconstruction Methods[END_REF] (also derived in [START_REF] Kudo | Derivation and implementation of a conebeam reconstruction algorithm for nonplanar orbits[END_REF]). The corresponding intermediate function G 1 will be denoted G s . Moreover, it is shown in [START_REF] Clackdoyle | Overview of reconstruction algorithms for exact cone-beam tomography[END_REF] that the method of Tuy [START_REF] Tuy | An Inversion Formula for Cone-Beam Reconstruction[END_REF] is mathematically equivalent to Grangeat's method and fits the case τ = 0. Though similar in their mathematical treatment, the two methods are very different in practice. As already encountered in the inversion of the Radon transform, the ramp filter ε 1 has infinite support and evaluation of G s requires that the whole projection of the object is available. On the other hand, convolving with δ acts as a differentiation and is thus a local operation. The consequences of these differences will be studied in Chapter II-1.

Though of theoretical importance, these methods have not been widely used in tomographic systems. One of the reasons is that the final step -the outer integral in Equation 64 -is not a backprojection. Instead, it scans all planes passing through the point being reconstructed, find a source position inside that plane, and add the corresponding contribution to the reconstruction. So methods following this approach require that all projection data must be available before reconstructing. In a (filtered) backprojection, each projection can be treated as it is acquired. FBP-like cone-beam reconstruction methods have been proposed, for example in [START_REF] Defrise | A cone-beam reconstruction algorithm using shift-variant filtering and cone-beam backprojection[END_REF][START_REF] Katsevich | Theoretically exact filtered backprojection-type inversion algorithm for spiral CT[END_REF][START_REF] Kudo | Derivation and implementation of a conebeam reconstruction algorithm for nonplanar orbits[END_REF].

Much more could be said about tomographic reconstruction. Specifically, the question of dealing with incomplete data (limited angle, truncated data) has not been addressed. The reader is referred to [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF] for a mathematical treatment of these questions. Also, we did not mention helical trajectories. They are used in most of the diagnostic CT systems today.

T H E G R A N G E AT T H E O R E M A N D I T S P R O O F

The previous section introduced the intermediate functions of Equation 55 as they will play a important role in the contributions presented in Part II. In this section, we give the Grangeat theorem (central to Chapters II-1 and II-3) in its original form. The proof starts from the intermediate function G g of Equation 55, with τ = 1. The proof relies essentially on the change of variables from spherical coordinates to Cartesian coordinates of the detector, which we derive completely here and which will be used several times in the sequel.

We consider a cone-beam projection g λ acquired on a detector, denoted D. For β ∈ S 2 we let H β, s λ • β denote the plane orthogonal to β and containing s λ . Without loss of generality, the DCS origin is at the principal point 10 and the u-axis is parallel to D H β, s λ • β . The 2D DCS is extended to 10 The principal point is the orthogonal projection of the source position onto the detector. a 3D coordinate system by including the unit vector w . In the 3DCS, we use (θ, φ) to denote the spherical coordinates, so defined that α(θ, φ) = (sin φ, sin θ cos φ, -cos θ cos φ), as illustrated in Figure 8 right. We also define θ so that β = (0, sin θ, cos θ). Finally, f denotes the distance from the source to the detector. Then we have Theorem 7 (Grangeat's Theorem).

1 cos 2 θ ∂ ∂v +∞ -∞ g λ (u, v) f u 2 + v 2 + f 2 du v= f tan θ = ∂ ∂s Rµ β, s s= s λ • β (65)
In the original work of Grangeat, the theorem was expressed in the form of Equation 65. The computation of the LHS from the projection image g λ (u, v) follows the following four steps:

1. Weight the projection with the cosine of the incidence angle, Proof. The proof proceeds by proving that both the LHS and the RHS are equal to the intermediate function G g λ, β . From Equation 62, we already have

G g (λ, β) = +∞ -∞ δ s -s λ • β Rµ β, s ds. ( 66
)
The right-hand side (RHS) of Equation 66is exactly the RHS of Equation 65. We will now prove that

O u v w s λ β α θ φ (u, v)
G g λ, β = 1 cos 2 θ ∂ ∂v +∞ -∞ g λ (u, v) f u 2 + v 2 + f 2 du v= f tan θ ( 67
)
by changing the integration over the sphere S 2 in G g to the integration over the detector's (u, v) coordinates. The first step consists in parametrizing the sphere S 2 with spherical coordinates. In the 3DCS (O, u, v , w ) defined above, we have α = (sin φ, sin θ cos φ, -cos θ cos φ) (68)

d α = cos φ dφ dθ (69) (θ, φ) ∈ [-π, π[ × ]-π/2, π/2] ( 70
)
We do not re-write the integral for now since we are not interested in having it in terms of spherical coordinates but keep in mind the cos φ term that arises from the parametrization. Instead, we move forward to the (u, v) coordinates with the following change of variables:

v = f tan θ (71) u = v 2 + f 2 tan φ ( 72
)
The three following facts hold:

dθ dφ = f (u 2 + v 2 + f 2 ) v 2 + f 2 d u d v, ( 73
)
cos φ = v 2 + f 2 u 2 + v 2 + f 2 , ( 74
)
α(θ, φ) = u u + v v -s λ u 2 + v 2 + f 2 (75)
Equation 73 results from differentiating Equations 71 and 72. Equation 74 results from simple trigonometry. Equation 75 results from the choice of the coordinate system ( u, v , w ) (see Figures. 8 and9). Putting all this together, slightly simplifying and using the homogeneity of δ (of degree -2), we have :

G g λ, β = δ u u + v v -s λ u 2 + v 2 + f 2 • β g λ (u, v) f du dv (u 2 + v 2 + f 2 ) u 2 + v 2 + f 2 (76) = δ (u u + v v -s λ ) • β g λ (u, v) f u 2 + v 2 + f 2 du dv (77)
Note that the remaining weight in the integral is the usual cosine of the incidence angle. We denote g the cosine weighted projection (this notation will be used in the rest of the document):

gλ (u, v) = f u 2 + v 2 + f 2 g (λ, u, v). ( 78
)
and remark that u • β = 0 by definition of u, that v • β = cos θ, and since s λ = f w , s λ • β = f sin θ (see Figure 9). We then have:

G g λ, β = R 2 δ v cos θ -f sin θ gλ (u, v) du dv, ( 79
) = R 2 1 cos 2 θ δ v -f tan θ gλ (u, v) du dv, ( 80
) = 1 cos 2 θ ∂ ∂v gλ (u, v) du v= f tan θ , (81)
which completes the proof. I2

D E S C R I B I N G T H E G E O M E T R Y O F O N E P R O J E C T I O N
This chapter focusses on the cone-beam geometry and presents a parametric model for the accurate description of the geometry of one conic projection. By geometry of one projection, we mean position and orientation of the gantry, i. e. the system made of the X-ray source and the 2D flat detector. Mapping 3D space points to 2D detector's pixels relies on the geometry hence the importance of this geometric description. We make the following hypotheses, which correspond to the physical devices that this work is addressing:

• The X-ray source is infinitesimally small,

• The surface of the detector is a rectangle, with length l and width w (expressed in millimetres),

• The pixels are squares (the grid is orthogonal and the size of the pixels is the same in both directions).

• The detector is equipped with a detector coordinate system (DCS), denoted (O d , u, v ). O d is the centre of the detector, the u-and v-axes coincide with detector rows and columns respectively.

N I N E G E O M E T R I C P A R A M E T E R S

With the above assumptions, a complete description of the geometry of one projection requires nine geometric parameters. We adopt the parametrization of the Computer Vision community to model a pinhole camera. In that convention, extrinsic parameters are distinguished from intrinsic parameters. The former describe the position and the orientation of the device. The latter describe the relative position of the two components (the X-ray source and the detector). Given a 3D patient coordinate system (PCS) (O, x, y, z), the source position is described by three parameters (x λ , y λ , z λ ) (where λ ∈ Λ is a scalar which parametrizes the trajectory) and the orientation of the detector axes is determined by three Euler angles (η, θ, φ) (see Figure 11 and below). Note that (θ, φ) alone, if seen as spherical coordinates as in Section 1.7, define a unit vector normal to the detector plane. The angle η finally rotates the detector axes within this plane. We now place the detector in space so that in the 3D detector coordinate system (3DCS), the source has coordinates (u 0 , v 0 , f ). Figure 10 recapitulates all the parameters.

O

N E U L E R A N G L E S A N D R O TAT I O N M AT R I C E S
Defining the 3D orientation of the detector is a problem of convention. Intuitively, it requires 23 In the 3DCS (O d , u, v , w ), the source has coordinates (u 0 , v 0 , f ).

s λ =     x λ y λ z λ     O O d v 0 u 0 e z e x
three degrees of freedom: two parameters define the normal direction of the detector. The third parameter rotates the detector into the plane defined by the normal. In summary, orienting the detector is equivalent to applying a 3D rotation to the axes of the PCS. A 3D rotation can be decomposed into the product of three elemental rotations 11 . There are 24 possible 11 We call elemental rotation a rotation about one axis of a coordinate system.

conventions for this decomposition 1 . One usually refers to such a convention as Euler angles. The convention which is used in this work is the Tait-Bryan convention with the sequence zxy. It is described in Figure 11. Any 3D rotation, denoted R η,θ,φ can be decomposed as R η,0,0 R 0,θ0 R 0,0,φ where each rotation is about the z-, x-and y-axis of each intermediate rotated coordinate system 12 . In matrix notation, we have:

rotations are defined by their action on the coordinates rather than with their action on affine points.

R η,θ,φ =     cos η -sin η 0 sin η cos η 0 0 0 1         1 0 0 0 cos θ -sin θ 0 sin θ cos θ         cos φ 0 sin φ 0 1 0 -sin φ 0 cos φ     (82)
This parametrization describes the acquisition geometry of one projection. The modification of one intrinsic parameter u 0 , v 0 or f without changing the eight other parameters results in the modification of the 3D position of the detector, without affecting the source position, nor the orientation of the detector. The modification of one of the source parameters x λ , y λ or z λ without changing the eight other parameters results in the modification of the whole system (source and detector). Regarding the three Euler angles, there is an ambiguity in what should be the state of the system after one of these angles is modified (without changing the eight others). In the Computer Vision world, this would result in a displacement of the detector (both in orientation AND position). But the mechanical effect of a deforming C-arm that we want to model does not fit this Computer Graphics convention. When the detector is tilted, we mean to change the orientation of the detector but not its position. The ambiguity is: what point of the detector is the reference point for its 3D position?

2.2 T H E P R O J E C T I O N M AT R I X
The acquisition geometry is a key input to the tomographic reconstruction.

In CT software, the geometry may be defined by an end-user using geometric parameters as described above. This information is then encoded in the system as a mapping from 3D point in world (homogeneous) coordinates to the 2D (homogeneous) coordinates of the central projection of that point on the detector, e. g. the intersection of the detector and the line passing through the source and the point. It is usual to encode this mapping as a matrix, called the projection matrix. This section explains how the projection matrix is built from a set of parameters. The reader who is not acquainted to projective geometry will find the basics in Appendix 1.3. The projection matrix is a projective mapping from a 3D-real projective space to the 2D-real projective plane of a detector. The size of the matrix is 3 × 4 (points are expressed in homogeneous coordinates).

The intrinsic matrix

Let us, for now, make the following assumptions on the geometry:

• The source is at the origin O of the PCS.

• The detector lies in the plane z =f ( f > 0).

• The origin O d of the detector is at (0, 0,f ).

• The u-and v-axes are parallel to and have same direction as the x-and y-axes respectively.

With these assumptions, it easily seen (see Figure 12) that a 3D point (x, y, z) projects on a point (u, v) satisfying the relations

u = -f x z (83) v = -f y z , ( 84
)
which can be written in a matrix form:

    u v 1         -f 0 0 0 -f 0 0 0 1         x y z     , (85)
where means equality up to a non-zero multiplicative constant. If the detector is translated within the plane z =f by a vector (-u 0 , -v 0 ) (the principal point now has coordinates (u 0 , v 0 ) in the detector's coordinates), this relation becomes:

    u v 1     =     -f 0 u 0 0 -f v 0 0 0 1         x y z     (86)
The matrix in Eq 86 is called internal matrix or intrinsic matrix. It is usually denoted K in the rest of the manuscript. Let us now suppose that the whole system (source and detector) is translated by s = (x λ , y λ , z λ ) and rotated by R η,θ,φ . The projection matrix P is then obtained by right-multiplying K by the matrix corresponding to the rigid motion

P = K R η,θ,φ     1 0 0 -x λ 0 1 0 -y λ 0 0 1 -z λ     . ( 87
)
If unambiguous, the rotation matrix is simply denoted R instead of R η,θ,φ for ease of notation. In the computer vision literature, the decomposition of a projection matrix can take several forms, each of which conveys specific information.

P = K R I 3 | -s (88) = K R| -t (89) = A|P •4 , (90)
where P • j denotes the j -th column of the matrix P . Similarly, P T i • will be used to denote the i -th row of the matrix P . In Equation 88, I 3 denotes the identity matrix of R 3 and s the source position in the PCS. In Equation 89, t denotes the coordinates of the source in the oriented coordinate system (OCS). The matrix R is the orthogonal matrix which converts PCS coordinates to OCS coordinates

R =     -e o x - -e o y - -e o z -     , (91)
so that t = R s. In Equation 90, the matrix A plays a significant role in the sequel. In the cone-beam context of this thesis, the matrix A is non-singular. The matrix A is singular if and only if the source position is at a point at infinity (in the 3D-real projective space), which corresponds to a parallel beam geometry.

F R O M T H E

P R O J E C T I O N M AT R I X T O T H E G E O M E T R I C P A R A M E - T E R S
In the previous section, we described how to build up the projection matrix from the geometric description of the system. In this section, we do the reverse. Given a projection matrix, how are the geometric parameters derived? The main mathematical tool is the so-called QR-decomposition. But some non-uniqueness issues in this decomposition would add irrelevant complexity to the presentation. Thus we derive the parameters manually.

It is easily seen that P •4 = -A s, hence the immediate identification of x λ , y λ and z λ :

    x λ y λ z λ     = -A -1 P •4 . (92)
Note that A is a 3 × 3 non-singular matrix whose inversion is not problematic in terms of computational load.

We now focus on the intrinsic parameters u 0 , v 0 and f . The rows of the matrix R are e o

x , e o y and e o z respectively. From the definition of A, we have

A =     -f 0 u 0 0 -f v 0 0 0 1         -e o x - -e o y - -e o z -     (93) =     --f e o x + u 0 e o z - --f e o y + v 0 e o z - - e o z -     . (94)
Using the orthogonality of the rows of R, we obtain

u 0 = A T 1• A 3• (95) v 0 = A T 2• A 3• ( 96
)
Recalling that we assumed f > 0, we also get

f = A T 1• A 1• -u 2 0 = A T 2• A 2• -v 2 0 ( 97
)
At this stage, the matrix K is fully determined. It follows that R = K -1 A. The last step is the derivation of the three Euler angles from the rotation matrix R. According to the zxy Tait-Bryan convention given in Equation 82, R is decomposed as

R =     cos η cos φ -sin η sin φ sin θ -sin η cos θ cos η sin φ + sin η cos φ sin θ sin η cos φ + cos η sin φ sin θ cos η cos θ sin η sin φ -cos η cos φ sin θ -sin φ cos θ sin θ cos φ cos θ     . ( 98
)
By choosing θ ∈ [-π/2; π/2] and the other two angles in [-π, π[, we have

θ = arcsin R 32 (99) φ = arctan2 (R 33 , R 31 ) (100) η = arctan2 (R 22 , -R 12 ) , (101)
where the function arctan2(x, y) returns the polar angle (in [-π, π[) of the point (x, y).

I3 D ATA C O N S I S T E N C Y C O N D I T I O N S : A R E V I E W

In this section, we give a general overview of the literature on data consistency conditions. From a very general perspective, let us consider a linear operator M : E -→ F which model the acquisition system (the measures).

If f ∈ E , the collected measures are g = M f . We address the following question: given g ∈ F , does g belong to the image (the range) of M ? DCC are conditions that the measured data must fulfill to belong to the range of the forward operator. Heuristically, if a function g does not belong to the image of M , the reconstruction process is doomed to failure, i. e. the reconstructed object function will be degraded by artefacts. Artefacts in image reconstruction from X-ray projections stem from physical effects (beam hardening, scattering...) which usually break the consistency of the measured data. A general roadmap for DCC-based artefact reduction is three-step:

1. Design a parametric model of the phenomenon (beam hardening, scatter, motion, misalignment...), 2. build a cost function which depends on the model parameters and quantifies the inconsistency of the data and, 3. minimize the cost function with respect to the model parameters.

T H

E R A N G E O F T H E R A D O N T R A N S F O R M
The Helgason Ludwig characterization is the most important. It provides a set of conditions which are necessary and sufficient. In short, these conditions give a complete answer to the problem of characterizing the range of the operator. After stating the main theorem and giving a partial proof, we give a geometrical interpretation of these conditions, which was proposed in [START_REF] Clarkson | Symmetry properties of an imaging system and consistency conditions in image space[END_REF]. I spend some space here to describe this approach for several reasons: first, I find the approach mathematically appealing and quite elegant. Second, this could provide a systematic roadmap for the derivation of geometrically driven DCC. I spent sporadic time on that topic over my PhD period and unfortunately was not able to elaborate such roadmap. It remains in the background of my to-do list.

The Helgason-Ludwig Consistency Conditions (HLCC)

The characterization of the range of the Radon transform has been given a complete solution in the mid 1960's by Ludwig [START_REF] Ludwig | The radon transform on euclidean space[END_REF] and Helgason [START_REF] Helgason | The Radon transform on Euclidean spaces, compact two-point homogeneous spaces and Grassmann manifolds[END_REF]. This characterization is now well known under the name of Helgason-Ludwig Consistency Conditions (HLCC).

Theorem 8 (Helgason-Ludwig Consistency Conditions). A function p is the Radon transform of a function µ ∈ S (R n ) if and only if:

1. p ∈ S S n-1 × R , 2. p is even in s and θ: ∀( θ, s) ∈ S n-1 × R, p -θ, -s = p θ, s , 3. For all integers m ≥ 0,

J m θ = R p θ, s s m ds ( 102
)
is a homogeneous polynomial 13 in the coordinates of θ of degree at 13 A polynomial

P (X , Y , Z ) is said homogeneous of degree m if it is a linear combination of monomials X i Y j Z k such that i + j + k = m. most m.
Proof. Most of the proof is omitted. Especially the sufficiency of conditions 1., 2. and 3.. We will not spend too much time on the first condition, which is technical. Our general assumption that the functions we are dealing with are in the Schwartz spaces is enough to guarantee the validity of most of the theorems covered in this work. The necessity of the second condition simply states that the Radon coordinates θ, s and -θ, -s define the same hyperplane: H β,s = H -β,-s . Finally, only the necessity of condition 3. remains. We simply substitute the definition of p θ, s in the definition of J m and apply the (unit-jacobian) change of variables x = s θ + y.

J m θ = R θ ⊥ µ s θ + y s m d y ds (103) = R n µ ( x) x • θ m d x (104)
The RHS of the last line is clearly 14 a homogeneous polynomial in the coor- 14 In 2D at least:

write θ = (cos φ, sin φ) and expand the x • θ m term.
dinates of θ.

In dimension two, the condition 3. is called in projection form. An expression will be said in projection form when its computation only involves one single projection. The condition 3. in projection form is equivalent to the scalar-product form:

∀0 ≤ m < |k|, 2π 0 R p θ φ , s s m e i kφ ds dφ = 0, ( 105
)
where θ φ = (-sin φ, cos φ) 15 . This condition is not in projection form since 15 This result generalizes to higher dimensions with spherical harmonics. See [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF] for instance.

all projections p( θ φ , •) are incorporated.

In the applications, the projection form will be preferred because the consistency of one projection can be checked as soon as it is acquired. On the other hand, the scalar-product form requires that all the projections are available to check their global consistency (see [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF] for more details and the section below). The reader is referred to the paper of Ludwig [START_REF] Ludwig | The radon transform on euclidean space[END_REF] for the original, complete, and general proof of Theorem 8. The proof of the sufficiency is quite involved and omitted here.

We now very briefly illustrate the polynomial nature of J m numerically. We simulated the 2D Radon transform of the 2D Shepp-Logan phantom over 180 degrees. The phantom and the corresponding sinogram are presented in Figure 13. We computed the functions J m (φ) for m = 0, 1, 2, 3. Results are reported in Figure 14.

A geometrical interpretation of the HLCC

This section is independent from the rest of the manuscript. It can be skipped without hindering comprehension of the following sections. It is directly inspired by [START_REF] Clarkson | Symmetry properties of an imaging system and consistency conditions in image space[END_REF]. In this paper, the authors present a geometrical approach for deriving DCC.

For ease of notation, we set E = S (R 2 ) and F = S (S 1 × R). E is called the object space and F is called the data space or Radon space. Elements of E are denoted µ (the object functions) and elements of F are denoted p (the projection functions). Both spaces are equipped with their standard inner-product

〈µ 1 , µ 2 〉 E = R 2 µ 1 ( x) µ 2 ( x) d x, ( 106
)
〈p 1 , p 2 〉 F = 2π 0 R p 1 φ, s p 2 φ, s ds dφ. ( 107
)
If θ φ = (-sin φ, cos φ), we consider the 2D Radon transform R :

E → F µ → Rµ (
n = 0 n = 1 n = 2 n = 3
Figure 14: The Helgason-Ludwig consistency conditions. For each m = 0, 1, 2, 3, the corresponding figure plots the moments J m (φ). The trigonometric nature of J m appears for each m. Note though that for the case m = 0, whether J 0 is a constant or not very much depends on the y-scale of the graphic.

such that [START_REF] Hertle | Continuity of the radon transform and its inverse on Euclidean space[END_REF] subspace of F . We let (imR) ⊥ be the orthogonal supplement of imR

Rµ θ φ , s = R δ s -θ φ • x µ ( x) d x. ( 109
) T H E R A N G E O F R A N D T H E N U L L S P A C E O F R * The range of R is denoted imR: imR = Rµ, µ ∈ E . It is a closed
(imR) ⊥ = p ∈ F such that 〈p, Rµ〉 F = 0 ∀µ ∈ E (110)
and null R * = p ∈ F such that R * p = 0 be the null space of the adjoint operator R * (null R * is a subspace of F ). We have the following Proposition 9. Let pi i ∈I be a basis of the null-space null R * . Then, p ∈ imR if and only if:

〈 pi , p〉 F = 0 ∀i ∈ I (111)
Before the proof, we have two lemmas. First,

Lemma 1. p ∈ imR ⇐⇒ 〈p, p〉 F = 0 ∀ p ∈ (imR) ⊥ (112)
Proof. The necessity of the condition is exactly the definition of (imR) ⊥ .

Suppose now that 〈p, p〉 F = 0 for all p ∈ (imR) ⊥ . Since F = imR ⊥ ⊕ (imR) ⊥ , one can decompose p as p R +p (imR) ⊥ with p R ∈ imR and p (imR) ⊥ ∈ (imR) ⊥ . The condition hence writes 〈p (imR) ⊥ , p〉 F = 0 for all p ∈ (imR) ⊥ which is equivalent to p (imR) ⊥ = 0. Hence p = p R and the proof is complete.

The following easy lemma relates (imR) ⊥ to the adjoint of R.

Lemma 2.

(imR

) ⊥ = null R * . ( 113
)
Proof. By definition of (imR) T

H E G R O U P O F D I R E C T S I M I L A R I T I E S S

We now introduce the general framework of [START_REF] Clarkson | Symmetry properties of an imaging system and consistency conditions in image space[END_REF]. We consider the set S of direct similarities of the plane R

s -1 = (1/r , -θ, -r -1 R -θ b). ( 115
)
Equipped with this internal law, S is a group.

G R O U P R E P R E S E N TAT I O N A N D S Y M M E T R Y G R O U P
A representation of S on E (or F or any vector space) is a linear group action of S on E , i. e. a map

S × E → E (s, µ) → s • µ , ( 116
) satisfying 1. ∀µ ∈ E , I • µ = µ 2. ∀(s, s ) ∈ S × S, s • s • µ = ss • µ 3. for λ ∈ R and µ, µ ∈ E , s ∈ S, s • (µ + λµ ) = s • µ + λs • µ .
For (s, µ) ∈ S × E , we define the representation of S on E by

s • µ ( x) = µ s -1 ( x) (117)
where s -1 ( x) is defined by Equations 114 and 115. It is easy to verify that Equation 117 defines a representation of S on E . In other words, for all s ∈ S, the map µ → s • µ, denoted π s , is an automorphism of E i. e. π s ∈ GL(E ). We denote π this representation. The representation π has the important property that 〈π s µ, π s µ 〉 E = a(s)〈µ, µ 〉 E for some non-zero a(s), for all s ∈ S and all µ, µ ∈ E . Such a representation is called a conformal representation.

Suppose now that we have two representations π and ρ of S on E and F respectively. In [START_REF] Clarkson | Symmetry properties of an imaging system and consistency conditions in image space[END_REF], S is called a symmetry group for R if R is an equivariant map, i. e. is such that the following diagram commutes

E F E F R π s ρ s R , (118)
which is equivalent to

R π s µ = ρ s Rµ , ∀s ∈ S, ∀µ ∈ E . ( 119
)
The representation π has already been defined. In order for S to be a symmetry group, for s = (r , θ, b) ∈ S, ρ s is defined by

ρ s (p) θ φ , s = p θ φ-θ , s -θ • b r . ( 120
)
It is also easily checked that 〈ρ s (p), ρ s (p )〉 F = r 〈p, p 〉 E so that ρ is also a conformal representation. So far, we have introduced two representations of the group S of direct similarities, which are compatible with the Radon transform R in the sense that Equation 119 is satisfied. Before moving forward to the derivation of the HLCC through this geometric machinery, it is necessary to emphasize an important property of (imR) ⊥ , the nullspace of R * (i. e. the orthogonal complement of imR):

Lemma 3. (imR) ⊥ is S-invariant: ∀s ∈ S, ∀ p ∈ (imR) ⊥ , ρ s ( p) ∈ (imR) ⊥ . ( 121
) Proof. Let s ∈ S and p ∈ (imR) ⊥ . Let µ ∈ E . One has 〈Rµ, ρ s ( p)〉 F = 〈ρ s ρ s -1 Rµ, ρ s ( p)〉 F (122) = r 〈ρ -1 s Rµ, p〉 F (123) = r 〈R π s -1 (µ) , p〉 F (124) = 0. ( 125
)
The second line is obtained with the conformality of ρ, the third one is Equation 119 and the last one is the hypothesis p ∈ (imR) ⊥ .

I N V A R I A N T S U B S P A C E A N D H L C C

Representation theory states that if a subspace is S-invariant, there must be a basis with functions of the form 16 :

16 This is the tricky part of the paper, which Iunfortunately -did not completely elucidate.

p (m,k) θ φ , s = s m e 2i πkφ , (126)
so that we recover the HLCC from Proposition 9.

The nice thing about this framework is that each part of the HLCC finds a nice geometric interpretation. The s m term in the basis functions is related to the scaling of the similarities. The exponential part is related to the rotation. And finally, the constraints on k and m in Equation 105 (0 ≤ m < |k|) result from a close analysis of translations. For other transform (and for 3D cone-beam in particular), the authors also derive DCC. We will mention them quickly at the end of Chapter II-1.

I N

D I V E R G E N T B E A M G E O M E T R Y
The HLCC are a complete answer to the consistency problem. Moreover, they are very useful in practice thanks to their projection form. When dealing with the divergent-beam geometry, it is a general rule that things get more complicated. As was already mentioned for the reconstruction methods, the first intuitive way of addressing a divergent-beam problem is to rearrange the divergent data into parallel-beam data. Of course, deriving DCC in divergent-beam geometry is no exception. We develop this approach in the first subsection, in dimension two. But the re-binning approach is not free from problems and fully divergent complete DCC were sought for. To the best of our knowledge, no complete DCC have been found so far for general acquisition geometry (without any restriction on the source trajectory). Even worse, no such set exist for the very useful acquisition geometry: an X-ray source orbiting around the object. Two situations can be found in the literature. First, necessary and sufficient DCC were derived at the price of quite severe restriction on the source trajectory. This will be addressed in the second subsection. Second, DCC were published for general trajectory, but with no information on their sufficiency. Some examples of such DCC will be presented in the third subsection.

The questions of sufficiency and independence of two sets of consistency conditions are of great interest though, in view of the applications. In the perspective of the general DCC-based artefact reduction method that we quickly sketched in the introduction of this chapter, a cost function is built by incorporating consistency conditions. The more independent information is incorporated in the cost function, the more parameters can be identified. If two sets of DCC are available, one would like to be able to determine whether the two sets of conditions are independent or not, in the sense that the information they each convey is independent from the other. If they are, exploiting both sets may increase the number of unknown parameters which can be recovered or increase the robustness of the method by enforcing extra independent constrains. If they are not, the evaluation of one set is useless. This problem is more precisely addressed in the third contribution of this work. We prove that two sets of DCC (the Grangeat Consistency Conditions [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF][START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF] and Fan-beam DCC presented in [START_REF] Levine | Consistency Conditions for Cone-Beam CT Data Acquired with a Straight-Line Source Trajectory[END_REF] and applied to the calibration problem in Chapter II-2), very different from each other in their theoretical derivations, are in fact closely related and not independent.

Rebinning fanbeam data on a circle to derive fanbeam DCC from HLCC

We recall the divergent beam transform D of an object function µ

Dµ λ, α = g λ, α = ∞ 0 µ s λ + t α dt . ( 127
)
The work by Finch and Solmon [START_REF] Finch | A characterization of the range of the divergent beam X-ray transform[END_REF] is an early attempt to rebin divergentbeam data to Radon data. In this work, they apply HLCC to the re-binned data to characterize the range of D. The essential practical limit of this work is the assumption on the source locus: projection data is assumed available for all source position on a sphere which contains the compact support of the object. This assumption guarantees that the re-binning operation will yield the complete Radon space data. Of course, it is extremely difficult to fullfill this requirement in a practical situation.

Closer to the practical constraints, we now restrict to fanbeam data and consider a circular trajectory of the source: λ ∈ [0, 2π[ denotes the angle that the source makes with the x-axis of the PCS, s λ = (cos λ, sin λ) (the radius of the trajectory is 1). Let g (λ, α) denote the divergent projection data of an object µ (whose support is assumed, without loss of generality, to be fully contained inside the trajectory of the source) acquired along that trajectory, with α ∈ S 1 . Remember that α can be parametrized in two common ways: 1. with φ, the angle made by α and the x-axis of the PCS. following re-binning formulas (see Section I-1)

g λ, φ = p φ, sin λ -φ , ( 128
)
g λ, γ = p λ -γ, sin γ , ( 129
)
where p = Rµ is the Radon transform of µ. Applying the corresponding changes of variables

s = sin λ -φ , ( 130
)
s = sin γ (131)
in the integral of the HLCC J m θ φ = R s m p θ φ , s ds, we get respectively

J m θ φ = g λ, φ cos(λ -φ) sin m (λ -φ) dλ, ( 132
)
J m θ φ = g φ + γ, γ cos γ sin m γ dγ. ( 133
)
Though appealing at first glance, the fanbeam to parallel-beam rebinning methods do not provide the nice feature of the HLCC that each newly acquired projection can be immediately checked for consistency 17 . 17 These remarks are extracted from [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF], II.B.

In the first case (Equation 132), the integration variable is the projection index, so that all projections are required. In the second case (Equation 133), the integration variable is the ray angle γ (which is convenient) but for each ray angle, the value to integrate should be picked up from a different projection (the projection g (φ + γ, •)).

Despite these drawbacks, the method was used in [START_REF] Patch | Moment conditions indirectly improve image quality[END_REF] to monitor defective detectors or in [START_REF] Kim | Data consistency-driven scatter kernel optimization for x-ray cone-beam CT[END_REF][START_REF] Patch | Moment conditions indirectly improve image quality[END_REF] to reduce the scattering artefacts in the reconstructed volumes. It was also used in [START_REF] Yu | Data Consistency Based Rigid Motion Artifact Reduction in Fan-Beam CT[END_REF][START_REF] Yu | Data consistency based translational motion reduction in fan-beam CT[END_REF] to estimate (and correct for) a rigid-motion of the patient in a fanbeam scan.

Necessary and sufficient consistency conditions for divergent-beam projections

The extension of HLCC in their projection form has been extended to the divergent beam geometries in two articles [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF][START_REF] Clackdoyle | Full data consistency conditions for cone-beam projections with sources on a plane[END_REF]. In the first one [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF], necessary and sufficient conditions are derived in the particular case of 2D fanbeam projections acquired from a source moving along a line. To some extent, the second paper [START_REF] Clackdoyle | Full data consistency conditions for cone-beam projections with sources on a plane[END_REF] is a generalization of the first one to 3D conebeam projections. In this extension, the X-ray source may vary within a plane instead of along a line. We start with the 2D case. The following 2D fanbeam DCC are crucial for the Chapters II-2 and II-3. The X-ray source moves along the line x = 0 in the (O, e x , e y ) PCS, so that s λ = (0, λ). We consider fanbeam projections of an object µ whose compact support is supposed to be strictly contained in the half-plane x > 0 (in particular, the support of µ does not intersect the source trajectory). Each ray is parametrized by the angle it makes with the x-axis (see Figure 17), so that the line integrals can be rewritten

Dµ λ, φ = +∞ 0 µ t cos φ, λ + t sin φ dt ∀φ ∈ - π 2 , π 2 . ( 134
)
The following theorem is published in [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF] Theorem 10 (Full fanbeam DCC for sources on a line). A function g (λ, φ) is the fanbeam projection of an object µ (g = Dµ) if and only if, for all m ∈ N the integral:

J m (λ) = π/2 -π/2 g λ, φ cos φ tan m φ dφ ( 135
)
is a polynomial in λ of degree at most m.

Proof. The proof of the sufficiency is based on a theorem by Edholm and is beyond the scope of this work. We focus on the proof of the necessity. We suppose g = Dµ for some object function µ and form the expression J m (λ).

s λ = (0, λ)• O α φ x • x y
J m (λ) = π/2 -π/2 g λ, φ cos φ tan m φ dφ, ( 136
) = π/2 -π/2 +∞ 0 µ t cos φ, λ + t sin φ dt tan m φ cos φ dφ. ( 137
) (138)
Polar coordinates (t , φ) are then changed to Cartesian coordinates according to (x, y) = (t cos φ, λ + t sin φ). The following easy facts hold:

t dt dφ = dx dy (139) tan φ = y -λ x . ( 140
)
Hence,

J m (λ) = µ x, y (y -λ) m x m dx dy x , ( 141
) = m k=0 c m,k λ k , (142)
with:

c m,k = (-1) k m k µ(x, y) y m-k x m+1 dx dy. ( 143
)
Since the support of µ is contained in the half-plane x > 0, there exists ε > 0 such that g (λ, φ) = 0 for |φ| > π/2 -ε. Hence, dividing by cos φ in Equation 135 does not introduce any singularity. The support condition, though, cannot be relaxed.

The conditions take the form of a polynomial in the projection index λ. The whole set of conditions are proved to be necessary and sufficient and this was the major contribution of the paper. The order-0 condition (J 0 is constant over the trajectory), on the other hand, had been known for decades and extensively exploited in many papers, along with other flavours of the same invariant integral conditions (see Section 3.2.3). The order-0 condition alone is of course not sufficient, but can be extended to more general trajectory, via pairwise DCC. Given a pair of projections acquired along any generic source trajectory, the order-0 can always be applied to the pair by considering the virtual linear trajectory that connects the two source positions. This remark is central in Chapters II-2 and II-3.

We now turn to the 3D cone-beam geometry. The same authors have extended their fanbeam result to the cone-beam acquisition geometry. In this 3D context, the source is constrained to lie in a plane (instead of a line) and the support of the object is assumed to not intersect this plane. The conditions take the same projection form, similar to HLCC. Let s λ = (x λ , y λ , 0) be the coordinates of the source in the PCS (the source trajectory is contained in the z = 0 plane). The object is supposed to be fully supported in the half-space z > 0. The projection data is the usual line integral

g λ, α = +∞ 0 µ s λ + t α dt (144)
With spherical coordinates (θ, φ), α = (cos φ sin θ, sin φ sin θ, cos θ) and for (X , Y ) ∈ R 2 , we define

J m (λ,U ,V ) = g λ, θ, φ U cos φ + V sin φ m tan m+1 θ cos θ dφ dθ (145)
Theorem 11 (Full cone-beam DCC for sources on a plane). A function g λ, θ, φ is the cone-beam projections of a compactly supported (in z > 0) object µ if and only if

• g has compact support for all λ,

• J m (λ,U ,V ) = K n (U ,V , -x λ U -y λ V ) for some homogeneous polyno- mial K m (X , Y , Z ) of degree m.
No proof is given here. For the proof, for the physical detector version of the theorem, and for linogram and planogram versions, the reader is referred to the original paper [START_REF] Clackdoyle | Full data consistency conditions for cone-beam projections with sources on a plane[END_REF]. Again, the support condition cannot be relaxed because the cos θ term in the denominator of the integrand of J m would introduce a singularity.

It should be noted that these conditions are truly cone-beam conditions in the sense that all the data of one projection is incorporated in the computation of J m . In this PhD work, I have not investigated the application

x y z • s λ x λ y λ O µ( x) φ θ
Figure 18: Cone-beam projections acquired from a source in a plane. The source is displaced in the plane z = 0. One ray is identified with the unit vector α with spherical coordinates (θ, φ) such that α = (cos φ sin θ, sin φ sin θ, cos θ).

of these DCC to the calibration problem addressed in Chapter II-2. The main reason is that they cannot handle truncated data (see condition 1 in Theorem 11), which is a typical situation in human body imaging. This limitation is also true for the linear fanbeam DCC of Theorem 10. But from the point of view of handling truncated data, fanbeam DCC extended to conebeam data allow to check the consistency of a limited, un-truncated subset of the fanbeam projection (typically those in the trans-axial direction).

Necessary only conditions in general acquisition geometries: the invariant integrals

The previous complete fanbeam DCC in Theorem 10 are an extension to any degrees m of some zero-order necessary conditions which have been known for much longer time. These conditions state the fact that the integral of a projection (weighted by some appropriate weight function) is independent of the projection index. The first version of the invariant integral version of the divergent-beam DCC dates back to the John's ultra-hyperbolic equation [START_REF] John | The ultrahyperbolic differential equation with four independent variables[END_REF] in a paper from 1938. In this section, the notation is borrowed from Patch's and Levine's paper [START_REF] Levine | Consistency Conditions for Cone-Beam CT Data Acquired with a Straight-Line Source Trajectory[END_REF][START_REF] Patch | Consistency conditions upon 3D CT data and the wave equation[END_REF] on the derivation of Consistency Conditions from John's equation [START_REF] John | The ultrahyperbolic differential equation with four independent variables[END_REF]. We let η = (η 1 , η 2 , 0) and ξ = (ξ 1 , ξ 2 , 1) denote coordinates of points in the planes Π 0 : z = 0 and Π 1 : z = 1 respectively. Points in Π 0 can be thought of as source positions and Π 1 as the detector plane. We denote g the weighted projection data with this parametrization of lines

g ξ, η = R µ ξ + t η -ξ dt . ( 146
)
We have the following relation with the usual divergent projection data

g ξ, η = 1 η -ξ g ξ, η -ξ η -ξ . ( 147
)
Note that the term 1/ ηξ is the cosine if the incidence angle. The projections g satisfy the John's ultra-hyperbolic equation 1818 See [START_REF] Patch | Consistency conditions upon 3D CT data and the wave equation[END_REF] for a proof.

∂ 2 ∂η 1 ∂ξ 2 - ∂ 2 ∂η 2 ∂ξ 1 g ξ, η = 0. ( 148
)
The original idea of using John's equation is to use the measured data as boundary conditions to infer non-measured data, in order to improve reconstruction quality. Instead, we are interested here in using John's equation to derive DCC. Taking the Fourier transform of Equation 148 with respect to η yields

κ ⊥ • ∇ ξ ĝ ξ, κ = 0, (149)
where κ denotes the Fourier dual variable of η and ĝ denotes the Fourier transform of g with respect to the η variable. Assuming now (without loss of generality) that the coordinates are such that ξ = (ξ, 0, 1) and κ ⊥ = (κ, 0, 0), the following consistency condition is derived from Equation 149:

R g ξ, η dη = R g ξ + t , η dη (150)
The last equation is a data consistency condition. The situation is depicted in Figure 19. Given two cone-beam projections with source positions on a line l , the integral of the projections weighted by the cosine of the incidence angle and restricted to a line parallel to l are equal. We now quickly establish the link between John's consistency condition of Equation 150 and the m = 0 case of the fanbeam DCC on a line of Theorem 10. In Figure 20, we represent the plane of interest of Figure 19 in the same configuration as in Figure 17. We then simply change the projection index ξ to the usual λ and turn the detector's coordinate η to the usual angular variable φ with: which is exactly the m = 0 case of the fanbeam DCC. In the rest of the manuscript, the acronym FBCC will refer to these particular zero-order fanbeam DCC. The same DCC have been published in various guises and derived through different routes.

η = ξ + d tan φ. Further noticing that ξ -η = (0, ξ) -(d , η) = d / cos φ, we obtain R g ξ, η dη = π 2 -π 2 g λ, φ d cos 2 φ dφ, (151)
• From a general formula established by Hamaker et al. [START_REF] Hamaker | The divergent beam x-ray transform[END_REF], Noo established a re-binning formula: the Hilbert projection equality (see Equation 10 in [START_REF] Noo | Image reconstruction from fan-beam projections on less than a short scan[END_REF] and Equation 15 in [START_REF] Clackdoyle | Tomographic Reconstruction in the 21st Century[END_REF]), from which an integral invariant DCC is easily derived.

• In [START_REF] Wei | Integral invariants for computed tomography[END_REF], the authors use a homogeneous function along with its group of symmetries to derive DCC in the same form of invariant integrals.

• Chen and Leng [START_REF] Chen | A new data consistency condition for fanbeam projection data[END_REF] obtained a similar equality, but their derivation was carried out essentially in the Fourier space.

In [START_REF] Tang | The mathematical equivalence of consistency conditions in the divergent-beam computed tomography[END_REF], all these invariant integral DCC were proved to be mathematically equivalent.

Grangeat-based DCC

The DCC presented in the previous section are essentially fanbeam consistency conditions. Even in a cone-beam geometry, the DCC of Theorem 10 are obtained by restricting the projection data to planes which contain a subset of source positions. In this section, we introduce Grangeat Consistency Conditions (GCC). They are based on the theorem of Grangeat (see Section 1.7 and [START_REF] Grangeat | Analyse d'un systeme d'imagerie 3d par reconstruction a partir de radiographies x en geometrie conique[END_REF][START_REF] Grangeat | Mathematical framework of cone beam 3D reconstruction via the first derivative of the radon transform[END_REF]) and were introduced in the context of jitter correction for cone-beam CT in [START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF]. The original idea was extended in [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF].

Let g λ, α be a set of cone-beam projections with source positions s λ , for λ ∈ Λ. Let β ∈ S 2 a unit vector and H β,s a plane orthogonal to β at signed distance s from the origin of the PCS. For a given projection, we define detector coordinates and spherical coordinates as described in Section 1.7 and recalled in Figure 21. We denote l the line of intersection of the plane H β,s with the detector. The Grangeat theorem states that In Section 1.7, Equation 67, we proved that the LHS of Equation 153 is equal to the intermediate function G g λ, β . The Grangeat-based consistency condition relies on the simple fact that if H β,s contains two source positions s λ 1 and s λ 2 , then the LHS must equate:

1 cos 2 θ ∂ ∂v +∞ -∞ g λ (u, v) f u 2 + v 2 + f 2
G g λ 1 , β = G g λ 2 , β ,
because the RHS do. The 1/ cos 2 θ can be dropped, which yields the following equivalent consistency condition:

Proposition 12 (Grangeat Consistency Conditions (GCC)). For all β ∈ S 2 , if s λ 2 • β = s λ 2 • β, then ∂ ∂v +∞ -∞ g λ 1 (u, v) f u 2 + v 2 + f 2 du v= f tan θ = ∂ ∂v +∞ -∞ g λ 2 (u, v) f u 2 + v 2 + f 2 du v= f tan θ (154)
The GCC differ from the invariant integrals presented in the previous section in that they are truly cone-beam. The differentiation step will require the projection data not only along the line l but also along another line l ε , close and parallel to l .

Another important remark is in order. The computation of the GCC is three steps: weight the projections with the cosine of the incidence angle, integrate on the u-direction and differentiate in the v-direction. It is necessary to emphasize that the direction of the derivative is orthogonal to the line of intersection of the plane H β,s with the detector. Since all the computation steps are independent of the particular epipolar geometry, things can be pre-computed once for each projection. This point will be discussed in the Chapter II-3.

Applications of DCC

DCC are applied to various artefact correction methods. In Chapter I-4, we present geometric calibration methods in general and consistency-based methods [START_REF] Basu | Feasibility of tomography with unknown view angles[END_REF][START_REF] Basu | Uniqueness of tomography with unknown view angles[END_REF][START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF][START_REF] Maass | Geometrical Jitter correction in computed tomography[END_REF] in particular. DCC are also used in beam-hardening correction [START_REF] S Abdurahman | Calibration Free Beam Hardening Correction Using Grangeat-Based Consistency Measure[END_REF][START_REF] S Abdurahman | Beam hardening correction for bi-material objects using Grangeat-based consistency measure[END_REF][START_REF] Mou | Beam hardening correction based on HL consistency in polychromatic transmission tomography -art. no. 691333[END_REF][START_REF] Tang | Data consistency condition-based beam-hardening correction[END_REF] and attenuation correction in PET [START_REF] Welch | Attenuation correction in PET using consistency information[END_REF] and SPECT [START_REF] Mennessier | Attenuation correction in SPECT using consistency conditions for the exponential ray transform[END_REF].

Another possible application of the DCC is the extrapolation of missing projection data by enforcing a consistency constrain on the extrapolated data set. This is proposed in [START_REF] Van Gompel | Elliptical extrapolation of truncated 2D CT projections using Helgason-Ludwig consistency conditions[END_REF] and [START_REF] Kudo | Sinogram recovery with the method of convex projections for limited-data reconstruction in computed tomography[END_REF] for instance, with the HLCC. Though, truncated data is a major issue regarding DCC. All the DCC presented in this chapter assume un-truncated data 19 . And in most clinical 19 In Chapter II-2 we explain that this constraint can be slightly relaxed in the case of invariant integrals like FBCC.

applications, projection data are truncated. An attempt in this direction is [START_REF] Clackdoyle | Data consistency conditions for truncated fanbeam and parallel projections[END_REF] but the extension of all the above artefact correction methods to truncated projections is still pending.

I4 G E O M E T R I C C A L I B R AT I O N O F T O M O G R A P H I C S Y S T E M

The sections 4.1 and 4.2 contain paragraphs extracted from our publication [START_REF] Lesaint | Calibration for Circular Cone-Beam CT Based on Consistency Conditions[END_REF], Section II.C.

A computer implementation of any of the reconstruction algorithms described in the previous section requires a precise mapping of the 3D coordinates of the patient coordinate system to 2D coordinates of the detector coordinate system, in order to back-project pixel values of the projections to the correct voxels in the reconstructed volume. This mapping usually takes the form of a projection matrix (see Section 2.2), for each acquired image. The geometric calibration of a tomographic system is the process which determines the projection matrices of each projection of a scan.

Much work has been done on the calibration of CT systems. We give a quick review of the methods and briefly summarize the relative importance of each parameter with respect to their impact on the reconstruction quality. The design of the system is the first answer to the misalignment problem: it may guarantee that the nominal geometric parameters of the system are kept -during a scan -within margins which do not impact the reconstruction quality. To account for mechanical shift overt time, it is necessary to periodically update the calibration, by performing new calibration procedures. These procedures may be sensor-based or imaging-based. A sensor-based calibration is completely out-of-scope of this thesis. It relies on the use of a Coordinate Measuring Machine (CMM) which physically senses various points of the CT system to determine their locations and derive the whole geometry. Of course, it is dependant on a proper calibration of the CMM itself. The imaging-based methods, far more popular in all CT devices, rely on the analysis of the images of the system itself. Imagingbased calibration methods fall into two broad categories: off-line and online methods, which are described in the following sections.

P R O J E C T I O N M AT R I C E S C A L I B R AT I O N V S . G E O M E T R I C P A R A M E T E R S C A L I B R AT I O N
The final product of the calibration process is the stack of projection matrices (one for each projection). A brute force method consists in estimating each of these matrices, i. e. 11 degrees of freedom 20 . 20 Remember that a projection matrix is a 3 × 4 matrix defined up to a scaling constant due to its projective nature. Hence 12-1=11 degrees of freedom.

We call this method brute force because it does not take advantage of certain geometric characteristics of the system which may simplify the process. For instance, when a numeric flat-panel detector is used, the pixel grid is known to be orthogonal with square pixels. This simple fact kills two degrees of freedom in the calibration problem but not two coefficients in the projection matrix. Another important example is that of a turning table micro-CT system like the one which is calibrated in Chapter II-2. The source and the detector are fixed during the scan of an object. Their position and orientation are independent of the projection. This reduces the number of parameters of the calibration problem from 9N (N is the number of projections) to 7 + N (see Chapter II-2 for details). Hence the need to use a parametric model of the geometry. This model was described in Section 2.1.

4.1 O F F -L I N E G E O M E T R I C C A L I B R AT I O N M E T H O D S
Off-line calibration techniques encompasses all the techniques which require the scan of a dedicated calibration phantom prior to the production scan. The outcome of such techniques has a limited lifetime, depending on the mechanical stability of the system. The calibration procedure is launched periodically (typically once a year for a conventional medical scanner). Such methods rely on a major assumption that the geometric configuration is reproducible over time: the geometric configuration at the time of the calibration scan is the geometric configuration at the time of a future production scan. Various methods have been proposed [START_REF] Cho | Accurate technique for complete geometric calibration of cone-beam computed tomography systems[END_REF][START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF][START_REF] Gullberg | Estimation of geometrical parameters and collimator evaluation for cone beam tomography[END_REF][START_REF] Mennessier | Direct determination of geometric alignment parameters for cone-beam scanners[END_REF][START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF][START_REF] Rizo | Geometric calibration method for multiple-head cone-beam SPECT system[END_REF][START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF]: they all require pre-scanning of a calibration phantom, usually made up of small radio-opaque markers whose relative positions are accurately known. The theoretical projections of the markers (which depend on the geometric parameters) are afterwards compared with their actual projections to derive -iteratively [START_REF] Gullberg | Estimation of geometrical parameters and collimator evaluation for cone beam tomography[END_REF][START_REF] Rizo | Geometric calibration method for multiple-head cone-beam SPECT system[END_REF] or analytically [START_REF] Cho | Accurate technique for complete geometric calibration of cone-beam computed tomography systems[END_REF][START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF][START_REF] Mennessier | Direct determination of geometric alignment parameters for cone-beam scanners[END_REF][START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF][START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF] -the calibration parameters. In [START_REF] Mennessier | Direct determination of geometric alignment parameters for cone-beam scanners[END_REF], they solve the complete calibration problem and analytically derive all nine parameters for each projection.

O N -L I N E G E O M E T R I C C A L I B R AT I O N M E T H O D S

The other group of techniques consists of on-line techniques. All methods in this category solve the calibration problem without a specific calibration scan of a calibration object. They only use the data from the projections of the imaged object. Beyond this common feature, this group encompasses substantially different techniques. In [START_REF] Kingston | Reliable automatic alignment of tomographic projection data by passive auto-focus[END_REF][START_REF] Kyriakou | Simultaneous misalignment correction for approximate circular conebeam computed tomography[END_REF][START_REF] Muders | Stable and Robust Geometric Self-Calibration for Cone-Beam CT Using Mutual Information[END_REF], they minimize a cost function, whose evaluation requires the reconstruction of the object from the current estimate of the geometric parameters. The metric is based on entropy in [START_REF] Kyriakou | Simultaneous misalignment correction for approximate circular conebeam computed tomography[END_REF], the L 2 -norm of the image gradient in [START_REF] Kingston | Reliable automatic alignment of tomographic projection data by passive auto-focus[END_REF] or the mutual information between re-projected image and projection data in [START_REF] Muders | Stable and Robust Geometric Self-Calibration for Cone-Beam CT Using Mutual Information[END_REF]. The limitation of such methods is the computational load, which may not fit clinical workflow (though [START_REF] Kingston | Reliable automatic alignment of tomographic projection data by passive auto-focus[END_REF] limits this drawback by only reconstructing a fraction of the volume). In [START_REF] Ouadah | Self-calibration of cone-beam CT geometry using 3D-2D image registration[END_REF], they use the 3D reconstruction of a planning CT and compute projection-specific geometric parameters by registering the actual projections with the re-projected CT image.

Other works in this category utilize the redundancy of the projection data (i. e. the DCC). In the 2D parallel beam case, Basu and Bresler [START_REF] Basu | Feasibility of tomography with unknown view angles[END_REF][START_REF] Basu | Uniqueness of tomography with unknown view angles[END_REF] solve uniquely and efficiently the problem of unknown projection angles and shifts with the Helgason-Ludwig DCC. Some works use the trivial "opposite-ray" condition [START_REF] Meng | Online Geometric Calibration of Cone-Beam Computed Tomography for Arbitrary Imaging Objects[END_REF][START_REF] Panetta | An optimizationbased method for geometrical calibration in cone-beam CT without dedicated phantoms[END_REF][START_REF] Patel | Self-calibration of a cone-beam micro-CT system[END_REF]. In [START_REF] Panetta | An optimizationbased method for geometrical calibration in cone-beam CT without dedicated phantoms[END_REF], this DCC, which normally only applies in the central plane (the plane of the trajectory), is extended to cone-beam projections of a particular class of symmetric 3D object functions and shows accurate calibration results when approximated in a central region of a generic object. More closely related to our work is a series of publications on epipolar consistency conditions [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF][START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF][START_REF] Maass | Geometrical Jitter correction in computed tomography[END_REF]. These DCC are based on the Grangeat theorem and relate the derivative of the 2D Radon transform of the projections to the derivative of the 3D Radon transform of the imaged object.

The comparison of previous works is not easy due to the parametrization which may differ with authors. Nevertheless, it is widely documented that the detector shift u 0 and the in-plane angle η are of crucial importance [START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF][START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF][START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF]. On the other hand, [START_REF] Muders | Stable and Robust Geometric Self-Calibration for Cone-Beam CT Using Mutual Information[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF] demonstrated that the two outof-plane angles (φ and θ) may be set to zero without affecting the image quality if their true values are kept below 2 • (which is a reasonable manufacturing accuracy requirement). Finally, miscalibration of the source-tocentre and source-to-detector distances does not introduce artefacts in the reconstructed volume and are therefore not calibrated. However, these two parameters affect the magnification of the reconstructed volume, which would not be acceptable in some cases, such as a metrology-oriented application.

Part II C O N T R I B U T I O N S

This part gathers the contributions of this PhD work. We present three contributions. We first introduce new zero-order cone-beam DCC for a circular trajectory (published in [START_REF] Lesaint | Two cone-beam consistency conditions for a circular trajectory[END_REF]) and present our contribution to their generalization to higher orders [START_REF] Clackdoyle | Data Consistency Conditions for Cone-Beam Projections on a Circular Trajectory[END_REF]. The second contribution is the calibration of a micro-CT system based on integral invariant fanbeam DCC (FBCC). This contribution was published in [START_REF] Lesaint | Calibration for Circular Cone-Beam CT Based on Consistency Conditions[END_REF]. The third contribution, not yet published, is an implementation of the integral invariant DCC in a physical flat-panel detector. This allows a theoretical comparison of FBCC and Grangeat-based DCC which seems to be new.

II1 C O N E -B E A M D C C F O R C I R C U L A R T R A J E C T O R Y
In the review of existing DCC (Chapter I-3, page 29), we mentioned essentially two cone-beam DCC: the Theorem 11 of Section 3.2.2 page 37 and the invariant integrals DCC (Section 3.2.3, page 40). The former addresses the special case of a planar trajectory with a planar detector parallel to the plane of the trajectory and the object placed in between these two planes (and in particular, not intersecting the trajectory plane). The latter are necessary conditions only. They can be used in a circular orbit cone-beam geometry but only in a pairwise fashion. This contribution introduces a new set of cone-beam consistency conditions for the very useful and common circular trajectory. The DCC are presented as they have been introduced (i. e. published): I was the first author of an abstract to the Fourth CT Meeting in Bamberg (Germany) in 2016 [START_REF] Lesaint | Two cone-beam consistency conditions for a circular trajectory[END_REF], which presented two modalities of a new zero-order DCC. Then, I contributed to the generalization of this result to higher degrees, which was published in the Signal Processing Letters in December 2016 [START_REF] Clackdoyle | Data Consistency Conditions for Cone-Beam Projections on a Circular Trajectory[END_REF]. I essentially contributed to the final version of the proof of the main result.

In this chapter, the X-ray source traverses a circle in the y = 0 plane, centred in O, with radius r , so that the source s λ has coordinates (-r sin λ, 0, r cos λ) in the PCS. Without loss of generality, we assume that the support of the object µ is fully contained in a sphere of radius r -ε for some ε > 0, so that the trajectory of the source never intersects the support of µ 21 . The detector is at distance f from the source and is aligned, i. e. the 21 Note that this assumption allows the object to intersect the plane of the trajectory.

principal ray passes through O (the origin of the PCS) and O d (the origin of the DCS) and the v-axis of the DCS coincides with the y-axis of the PCS (the axis of rotation). The Figure 22 illustrates the situation.

1.1 N E W Z E R O -O R D E R C O N E -B E A M D C C F O R A C I R C U L A R T R A J E C - T O R Y
The results presented in this section were published in J. Lesaint

G τ λ, β = S 2 ε τ α • β g λ, α d α, ( 155
)
where ε τ (s) = (1 -τ)δ (s) + τρ(s) 22 . We consider in this section the two ex- 22 The "functions" δ and ρ were defined in Chapter I-1 page 5. tremal cases: G 0 and G 1 , which we denoted G g and G s respectively. G s is the intermediate function proposed by Smith [START_REF] Smith | Image Reconstruction from Cone-Beam Projections: Necessary and Sufficient Conditions and Reconstruction Methods[END_REF]. G g was used by Grangeat [START_REF] Grangeat | Analyse d'un systeme d'imagerie 3d par reconstruction a partir de radiographies x en geometrie conique[END_REF]. The functions G g and G s are related to the 3D Radon transform of the object function via the relations 23 23 See Section 1.6.4.

G s λ, β = ρ * R β µ s λ • β (156) G g λ, β = -δ * R β µ s λ • β . ( 157
)
Consistency conditions can be inferred from any plane H β,s which contains at least two source positions s λ 1 and s λ 2 . In this case s λ • β = s λ • β so that (see Chapter I-3)

G s λ, β = G s λ , β , (158)
G g λ, β = G g λ , β . (159)
By setting β = e y , considering the particular trajectory plane H e y ,0 , and denoting Ḡ the corresponding functions of the projection index λ ( Ḡs (λ) = G s (λ, y) and Ḡg (λ) = G g (λ, y)), the pairwise consistency conditions propagate to the whole trajectory. We obtain the following cone-beam DCC: Proposition 13. The functions Ḡs (λ) and Ḡg (λ) are constant. They do not depend on the projection index.

We now turn to the practical implementation of these DCC and the respective expected behaviour of each filter ρ and δ . The computation of the Grangeat intermediate function with a physical detector has already been derived in Section 1.7, page 20, which in our particular case, becomes

Ḡg (λ) = ∂ ∂v +∞ -∞ f g λ (u, v) u 2 + v 2 + f 2 du v=0 . ( 160
)
As for the computation of the Smith intermediate function, a similar derivation leads to

Ḡs (λ) = +∞ -∞ ρ(v) +∞ -∞ f g λ (u, v) u 2 + v 2 + f 2 du dv. (161)
A few remarks are in order:

• The particular choice of the plane H y,0 allows a relatively simple evaluation of the function because the integration and convolution directions (u and v respectively) coincide with the pixel grid. On the contrary, a plane in general position will require some interpolation in the projection to compute the line integral.

• These DCC are truly cone-beam since their computation requires more projection data than the restriction of the projection to a single plane.

• Nevertheless, ρ and δ acts very differently with respect to the amount of cone-beam data incorporated in the functions Ḡ. The differentiation in Ḡg is local, whereas the convolution with ρ in Ḡs encompasses the entire projection because the support of the ramp filter is infinite. We will discuss this point below with the numerical simulations.

• The derivation of the link between G τ and the Radon transform (Equations 156 and 157 or Equation 62 and its derivation page 18) makes an essential use of the homogeneity of the distributions ρ and δ . Furthermore, any function with this homogeneity property must be a linear combination of ρ and δ .

• Any value of τ ∈ ]0, 1[ would have produced equally valid DCC. By mitigating ρ and δ , one would loose their specificity in terms of ability to detect inconsistencies (see below) and one would gain nothing since the DCC derived from such τ would not be independent of those derived from ρ and δ .

In order to numerically validate the Proposition 13 and to test the ability of the newly proposed DCC to detect inconsistencies of the projection data, we simulated the projections of a 3D Shepp-Logan phantom 24 with the Re- 24 The definition of the 3D Shepp-Logan phantom can be found in Kak and Slaney textbook [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF], on page 102.

construction Toolkit (RTK) [START_REF] Rit | The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK)[END_REF] software package 25 . We simulated 72 pro- jections, equally spaced over a 360 °angular range. The source-to-isocenter distance was set to r = 100 mm and the focal distance to f = 200 mm. The 3D Shepp-Logan phantom was scaled to semi-axes of (55.2, 73.6, 72) mm. The simulated projections were 1024 × 1024 pixels with a pixel size of 0.5 mm.

We focussed our experiments on two sources of inconsistencies: data truncation and motion of an object. In this experiment, we severely truncated 3 projections out of the 72 projections. Figure 23 illustrates the scenario. The projection 10 is truncated axially, but the projection data is complete around the central line. This situation is very common in a medical context where the entire patient body cannot fit the field-of-view of the scanner. The projection 40, on the other hand, is truncated on the right side. This kind of truncation impacts the central line. This corresponds to a large patient. Finally, an interior square of projection 55 was set to zero. Such truncation simulates a (massive) default in the detector or an occlusion. Truncating the projections in different positions with respect to the central line v = 0 allows for a better understanding of the respective response of the ramp-and derivative-based DCC. For each scenario, we evaluated both functions Ḡs (λ) and Ḡg (λ) at each projection index λ. Integrals were computed with a simple trapezoid quadrature rule. The derivative was implemented with a central finite difference. For the unmodified version of the projection data, i. e. for consistent data, we expected a constant signal for both functions. It would not be relevant to plot such signals with consistent data since the signals cannot be rigorously constant (due to numerical errors, essentially stemming from the quadrature and differential quadrature rules used in the evaluation of Ḡ). The constantness of the functions must be evaluated against an inconsistency which theoretically breaks the constant, as done in the following.

R E S U LT S

We analyse the results presented in Figure 25. We notice that the signals are not strictly constant. For example, between projections 11 and 39, the plots are slightly bumpy while the data are consistent. This point was mentioned above. We now examine the first experiment with truncated data (left column of Figure 25). The lateral truncation of projection 40 triggers a significant spike in both plots, which is higher in amplitude than the background noise. This enables the automatic detection of the inconsistency (by simply thresholding the signal for instance). The main difference between both functions Ḡg and Ḡs is revealed by observing projections 10 and 55. The simulated truncation in these two projections is not detected by Ḡg (no spike in the plot) while it is by Ḡs . The reason is that Ḡg only checks data which are close to the central line, while Ḡs scans the whole projection. In the first case, an inconsistency (truncation) which is away from the central line does not break the particular Grangeat condition. In the second case, the convolution with the ramp filter at v = 0 involves data over the support of the ramp, which is infinite.

The same kind of comment remains valid for the motion scenario. The Grangeat conditions detects motion inconsistency no sooner than one projection before the moving object intersects the central line, due to the central finite difference scheme to compute the derivative. The inconsistency introduces a sharp discontinuity in the consistency function. On the other hand, the DCC based on the ramp filter presents a smoother detection of the motion. Theoretically, the motion-induced inconsistency is detected as soon as it originates. In practice, the plot visibly deviates from the constant value not before projection index 40. The ramp filter decreases rapidly towards zero, so the amplitude of the response to the motion reduces as the moving object gets farther from the central line and it is mitigated with the background numerical noise. Depending on the problem which is being addressed, the use of the ramp filter or the derivative will be preferred. For instance, if the data are known to be axially truncated, the use of Ḡg seems more relevant since it will not be sensitive to structural inconsistency. On the contrary, if data are known to be un-truncated, the use of Ḡs will act as a finer comb as any inconsistency will be detected, regardless to where it originates.

One limitation is the background noise, which prevents detection of events of small amplitude (with respect to the two consistency metrics). To address this limitation, we applied a smoothing filter to the line integral signals of each projection. More precisely, the line integrals signals l (v k ) = g (u, v k ) du were smoothed with the formula

l (v k ) = 1 2 10 k+5 i =k-5 10 i l (v i ), (162)
for those k for which the formula makes sense. This filter consists in applying five times consecutively a three-point smoothing filter

l (v k ) = (l (v k-1 ) + 2l (v k ) + l (v k+1
))/4. The rationale for smoothing in the v-direction is that the quadrature rule is the source of numerical errors 26 .

[][-2cm]Smoothing directly the signal Ḡ

would make no sense, as it would mitigate data from different projections.

Both filters were applied to the smoothed signals and the results are presented in Figure 26. The effect of the smoothing kernel in the v-direction We have introduced two new DCC for a circular orbit cone-beam CT. These DCC are a rewording, in terms of DCC, of the intermediate functions introduced by Grangeat or Smith in their reconstruction algorithms. Though truly cone-beam, these DCC only apply in the central plane (i. e. the plane of the trajectory). We investigated the pros and cons of each of the DCC with respect to two sources of inconsistency -motion and truncation. The size of the support of each filter plays a major role and may serve differently the problem of interest. Finally, these DCC are only necessary conditions and certainly not sufficient. For instance, the Grangeat DCC cannot be sufficient as the projection data in a small neighbourhood of the central plane cannot determine uniquely the whole object. They can be thought of as zero-order DCC, which are extended to higher orders in the next section (though the sufficiency is still not addressed there).

1.2 N E W G E N E R A L C O N E -B E A M D C C F O R A C I R C U L A R T R A J E C T O R Y
The results presented in this section were published in R. Clackdoyle I was a coauthor of this paper and do not claim the original ideas therein. I essentially contributed to the proof of the main result of the paper and to some numerical validation. For self-completeness, a complete overview of the paper is given.

In this paper, we propose new cone-beam consistency conditions, which extend the Smith version of the DCC of Proposition 13 to higher orders. The acquisition geometry we consider here is exactly the same as in the previous section: a circular orbit of the source in the y = 0 plane and a detector undergoing the same trajectory, in opposite direction. The detector is perfectly aligned: in the PCS, the vectors ( e o

x , e o y , e o z ) of the OCS have coordinates (-sin λ, 0, cos λ), (0, 1, 0) and (cos λ, 0, sin λ) respectively. As will be described in the following, the proposed DCC involve singular integrals if the object intersects the plane of the trajectory. We will explain how this singularity can be handled by the use of distributions and smart numerical implementation.

We start from a modified version of the function Ḡs

Ḡm (λ) = S 2 g λ ( α)ρ m α • e o y α • e o x m d α, ( 163
)
where ρ m is a generalized version of the ramp filter ρ, which we define by its Fourier transform 27 27 Again, this definition must be understood in the sense of distributions.

ρ m (s) = R (-2i π) m+2 2(m + 2)! |σ|σ m e 2i πσs dσ. ( 164
)
Up to the (-π 2 ) multiplicative constant, ρ 0 is the usual ramp filter. A simple interpretation of ρ m can be given in terms of derivatives. The ramp filter ρ 0 is the successive application of a Hilbert transform (Fourier transform:

-i sign(σ)) and one single derivative (Fourier transform: 2i πσ). In a way, ρ m continues the process by further differentiating the relation n + 1 times (see below for the formal derivations and for a rigorous way of handling potential singularity which may arise). We already emphasized the importance of the homogeneity property of ρ. The distribution ρ m is also homogeneous, with degree -(m+2). This is easily seen by performing the change of variables σ = aσ in the expression ρ m (as), for a > 0.

The main result of the paper is two-fold. First, we prove that Ḡm is constrained to be a polynomial (the DCC) and second, we give the details for a numerical implementation of Ḡm in the detector. Proposition 14 (Cone-beam DCC). If g is the projection data of some object µ, then Ḡm (λ) is a homogeneous polynomial of degree m in (cos λ, sin λ), i. e.

Ḡm (λ) = m k=0 a m,k sin k λ cos m-k λ, ( 165
)
for some coefficients a m,k , k = 0, ..., m to be determined.

Proof. Substituting g λ with its definition, changing from polar coordinates (t , α) centred in s λ to Cartesian coordinates x centred in O (i. e. x = s λ + t α with t 2 dt d α = d x) and using the homogeneity of ρ m , we have

Ḡm (λ) = S 2 g λ ( α)ρ m α • e o y α • e o x m d α (166) = S 2 +∞ 0 µ s λ + t α ρ m α • e o y α • e o x m d α dt (167) = R 3 µ( x)ρ m x -s λ x -s λ • e o y x -s λ x -s λ • e o x m d x x -s λ 2 (168) = R 3 µ( x)ρ m ( x -s λ ) • e o y ( x -s λ ) • e o x m d x. ( 169
)
It is clear that

s λ • e o y = 0 = s λ • e o x , so that Ḡm (λ) = R 3 µ( x)ρ m x • e o y x • e o x m d x (170) = m k=0 m k (-1) k R 3 µ( x)ρ m ( x • e o y )x k z m-k d x sin k λ cos m-k λ. ( 171
)
The last expression is the desired polynomial form of Ḡm with

a m,k = (-1) k m k R 3 µ ( x) ρ m y x k z m-k d x. ( 172
)
We now turn to the implementation of Ḡm in the physical detector, i. e. we change from the variable α to the (u, v) coordinates:

α = (u e o x + v e o y - f e o z )/ u 2 + v 2 + f 2 .
We already did that change several times (see page 20), we just remind here the general formula, for any generic function b and with the usual notation

S 2 b( α) d α = b(u, v) f (u 2 + v 2 + f 2 ) 3/2 du dv. ( 173
)
Applied to the function Ḡm , we obtain

Ḡm (λ) = ρ m u e o x + v e o y -f e o z u 2 + v 2 + f 2 • e o y u e o x + v e o y -f e o z u 2 + v 2 + f 2 • e o x m g λ (u, v) f du dv (u 2 + v 2 + f 2 ) 3/2 , ( 174
) = gλ (u, v)ρ m (v)u m du dv. ( 175
)
where again the homogeneity of ρ m plays a major role. Very similar to the previous zero-order case, the computation consists in a pre-weighting of the projections, an integration in the u-direction followed by a filtering in the v-direction.

F O C U S O N T H E F I LT E R ρ m
All the derivations below are formal. They are legal when the object's support does not intersect the trajectory plane. If it does, we will explain how they can be justified. We start from the Hilbert kernel:

π -i sign(σ) e 2i πσs dσ = 1 s , ( 176
)
and iteratively differentiate both sides of the equality (usual differentiation in the RHS, multiplication by 2i πσ in the Fourier space LHS). We obtain:

-1 2 (2i π) 2 σsign(σ) e 2i πσs dσ = - 1 s 2 , ( 177
) -1 2 (2i π) 3 σ 2 sign(σ) e 2i πσs dσ = 2 s 3 , ( 178
) -1 2 (2i π) 4 σ 3 sign(σ) e 2i πσs dσ = - 6 s 4 , (179)
. . .

-1 2 (2i π) m+2 σ m+1 sign(σ) e 2i πσs dσ = (-1) m+1 (m + 1)! s m+2 . ( (180)
) 181
All these computations can be summarized in the following: ρ m (s) = 1/s m+2 . Of course, this is misleadingly simple since the Hilbert relation of Equation 176 is only valid for |s| > 0.

T H E O B J E C T D O E S N O T I N T E R S E C T T H E T R A J E C T O R Y P L A N E

We assume that µ( x) = 0 for all x = (x, y, z) such that |y| < ε, with ε > 0. In this case, g λ (u, v) = 0 for all |v| < ε. The integral in the definition of Ḡm (λ) can be decomposed as

Ḡm (λ) = |v|<ε gλ (u, v)u m du ρ m (v) dv+ |v|≥ε gλ (u, v)u m du ρ m (v) dv. ( 182
)
The first term in the RHS vanishes. In the second term, ρ m (v) can be substituted by 1/v m+2 , thanks to the derivations above, which become valid on |v| ≥ ε. We finally obtain the singularity-free DCC Proposition 15.

Ḡm (λ) = gλ (u, v) u m v m+2 du dv ( 183
)
is a polynomial in (cos λ, sin λ) of degree m. The integral in Equation 183 is easily implemented with two consecutive trapezoid quadratures.

I implemented this limited version of the DCC on 50 projections of a Shepp-Logan phantom, vertically offset so that the support of the phantom does not intersect the y = 0 plane. The projection with index 20 was axially truncated and half of the projection with index 40 was set to zero. The results 28 are reported in Figure 27. The trigonometric polynomial nature of 28 These simplified version of the DCC was not presented in the reference paper [START_REF] Clackdoyle | Data Consistency Conditions for Cone-Beam Projections on a Circular Trajectory[END_REF].

Ḡm is very well illustrated on the plots. The inconsistency introduced in projections 20 and 40 are clearly materialized as outliers. Though, the amplitude of the deviation from the consistent Ḡm is attenuated with higher m.

T H E O B J E C T D O E S I N T E R S E C T T H E T R A J E C T O R Y P L A N E

In the case where the support of the object intersects the y = 0 plane, the derivation above are not valid any more (because Equation 176 is not). The definition of Ḡm should be understood in the sense of distributions. In this framework, ρ m is the finite part fp 1/s m+2 distribution 29 . Its implementation 29 See Appendix1. 1. is done in the Fourier space by band-limiting the distribution ρ m and infer samples of ρ m in the direct space, before calculating the convolution in zero (see [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF], Section 3.3.3 for an implementation of the ramp filter ρ = ρ 0 , which can be generalized to ρ m ). The first naive implementation of the DCC was very unstable. The simulations presented in the paper were done by heavily regularizing the filter ρ m (by apodizing with a triangular window instead of a simple window function), and applying the same smoothing filter as the one described in Equation 162.

The simulations are illustrated in Figure 28. The theoretical trigonometric behaviour of Ḡm (λ) is graphically confirmed. However, this has been obtained at the price of heavy regularization and smoothing. This has two consequences. First, the singularity in the y = 0 plane seems difficult to handle properly. Second, the sensitivity of the DCC to an inconsistency will be accordingly lowered. And the previous experiment has shown that massive truncation is hardly detected in higher order DCC (without regularization nor smoothing). This instability limits the practical use of these DCC and explains why we did not push investigations further with the DCC.

A F U T U R E R E S E A R C H D I R E C T I O N

The DCC of Proposition 14 combined with its detector form of Equation 183 are in projection form. But they can easily and equivalently be rewritten as a scalar product. Indeed, being a homogeneous polynomial of degree m is equivalent to being orthogonal to the Fourier basis functions e i l λ with |l | > m, which we can write On the other hand, in Chapter I-3 (page 29), we described a geometrical framework inspired by a research article from 1997 [START_REF] Clarkson | Symmetry properties of an imaging system and consistency conditions in image space[END_REF], which used the symmetry group of the forward operator to derive consistency conditions. The paper applies this principle to several forward operators (the case of the 2D Radon transform was presented in Chapter I-3), among which the circular orbit cone-beam transform. In the paper, the axis of rotation is (Oz) and the elements of the unit sphere are parametrized by longitude and colatitude where l , q, m are integer indices satisfying 1 ≤ q ≤ m and l + m ≥ 0. Whether the consistency given in Equations 185 and 187 are related to each other remains to be investigated. Of course, in Equation 185, we have m = q hence a much smaller number of equations/conditions.

II2 C O N S I S T E N C Y -B A S E D C A L I B R AT I O N O F A C O N E -B E A M M I C R O -C T

This section reproduces the full version of our published paper J. Lesaint The Section 2.3.3 is the review of existing calibration methods. It is completely redundant with Chapter I-4, so it can be skipped with no arm for further understanding.

Sometimes, notation of the paper may conflict with notation previously established in this work. Nevertheless, the reading of this chapter is completely independent of the rest of the dissertation and entirely selfcontained.

The final section was not in the original paper. It is an addendum which further investigates the valley-shaped cost function in the plane of the parameters θ and v 0 .

A B S T R A C T

In cone-beam computed tomography (CT), imprecise knowledge of the acquisition geometry can severely impact the quality of the reconstructed image. This work investigates geometric calibration using data consistency conditions (DCCs). Unlike the usual marker-based off-line methods, the proposed method does not require the extra-scan of a calibration phantom. It is based on the minimization of a cost function, which measures the inconsistency between pairs of projections. The method has been applied to both simulated and real data. The latter were acquired from a micro-CT system with circular trajectory, for which the problem reduces to identifying global misalignments of the system. When compared to uncorrected reconstruction, the method significantly improved the image quality. When compared to marker-based calibration method, the image quality was similar but no calibration scan was required. Finally the method can handle axially truncated data. Axial truncation is very common in the medical context but often considered intractable for DCC-based methods. We also demonstrate DCC calibration from real data with axial truncation.

I N T R O D U C T I O N

In cone-beam computed tomography (CBCT), a 3D image is reconstructed from a set of 2D projections acquired from a point-like X-ray source. Poor 65 reconstructed image quality can arise due to many possible causes. One of these causes is an inaccurate calibration of the system. Calibration is the process through which the geometry of acquisition of the projections is accurately determined. By geometry of acquisition, we mean the position and orientation of the detector and the position of the X-ray source in a fixed reference frame.

Calibration of a CBCT system has been studied for a long time. Early works used the scan of a dedicated phantom to estimate the geometric parameters. These are known as off-line methods, to emphasize the need for a preliminary scan. The calibration scan provides accurate geometry information as long as the misalignments that were estimated are reproducible over time. In the extreme case where mechanical flexibility of the system makes the reproducibility assumption false, these methods become invalid. More recently, on-line (or self-) calibration methods have been developed. For each acquisition, and before reconstruction, the calibration is computed directly from the projections.

This work presents an on-line method based on the minimization of a cost function, that quantifies the inconsistency of the set of measured projections. The data consistency conditions (DCCs) that are incorporated in the cost function have been described in many different works (see Sec. 2.3.3). They are essentially fanbeam consistency conditions for a linear trajectory. They have been adapted to a circular trajectory CBCT system by re-sampling each pair of projections into a virtual detector parallel to the line connecting pairs of source positions. This idea was already proposed in [START_REF] Levine | Consistency Conditions for Cone-Beam CT Data Acquired with a Straight-Line Source Trajectory[END_REF] but, to the best of our knowledge, never implemented or applied to any CT reconstruction problem. Our work is very similar in its geometric approach to other recent works [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF][START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF] but differs fundamentally in the DCCs which are used. The proposed method was applied to simulated and real data, and compared with a robust off-line method.

N O TAT I O N A N D B A C K G R O U N D

The micro-CT system to which the calibration method was applied is made of a 2D flat detector and a fixed X-ray source. A turntable placed between the two allowed a full 360 • rotation of the object (see Figure29), so the acquisition geometry was equivalent to a circular trajectory of the source and detector. We will describe the geometry in detail, together with the calibration parameters that we are trying to estimate.

Description of the geometry

We use the same geometric description as that given in [START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF]. The detector cells are perfect squares (same width and height, arranged on a Cartesian grid). Let (x, y, z) be a fixed reference frame, defined as follows: the y-axis is the axis of rotation of the turntable. The origin is set so that the source lies in the y = 0 plane. The z-axis contains the source at rotation angle 0 and points in the direction of the source. The x-axis is defined so that (x, y, z) is a right-handed coordinate system. The flat panel detector is equipped with a direct (u, v) coordinate system whose origin is the center of the detector and whose axes coincide with the pixel rows and columns respectively. The geometry of one projection can be unambiguously described with eight parameters (see Figure 30):

• The rotation angle λ, taken from the z-axis,

• The radius of the source trajectory R,

• Three orientation angles of the detector (φ, θ, η),

• The source to detector distance D (or focal distance),

• The coordinates (u 0 , v 0 ) of the principal point (orthogonal projection of the source onto the detector plane).

With this parametrization, for λ ∈ [0; 2π[, the source position is given by s λ = (R sin λ, 0, R cos λ). The orientation of the detector is described with three Euler angles η, θ and φ (called yaw, pitch, roll respectively in [START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF] and skew, tilt, slant in [START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF]) applied in this order (respective axes of rotations are illustrated in Figure 30). The normal to the detector is defined with two outof-plane angles θ and φ about the u-and v-axes respectively. The in-plane rotation (about the focal axis) is given by η. The circular geometry thus consists of 8 degrees-of-freedom, unless the relative position and orientation of the source and detector can vary across projections. In our micro-CT system, the source and the detector are fixed, so the only projection-specific parameter is the rotation angle. The other seven parameters remain constant through the acquisition cycle. We call these parameters global misalignment parameters or global geometric parameters and refer to the corresponding geometry as true geometry. The system is perfectly aligned when (1) the principal axis (orthogonal to the detector plane and passing through the source) contains both the world origin and the detector origin and (2) the u and v axes of the detector are parallel to the x and y axes of the world frame at rotation angle λ = 0. In terms of the geometric parameters, these two conditions are equivalent to θ = φ = η = u 0 = v 0 = 0. We refer to the corresponding geometry as nominal geometry.

The X-ray line-integral model

If f ( x) = f (x, y, z) denotes the object density function, the projection g (λ, •) is defined by the usual line integral model:

g (λ, α) = ∞ 0 f ( s λ + t α)dt , ∀ α ∈ S 2 , ( 188
)
where S 2 denotes the unit-sphere of R 3 . The projection g (λ, •) vanishes for all α such that the line originating at s λ and directed by α does not intersect the support of f (see Figure31).

Review of existing calibration methods

Much work has been done on the calibration of CT systems. We give a quick review of the methods and briefly summarize the relative importance of each parameter with respect to their impact on the reconstruction quality.

If N denotes the number of acquired projections, the most general calibration problem consists of estimating -for each projection -the eleven Figure 31: Cone-beam geometry with circular trajectory. For a given scalar λ (typically in [0, 2π]), s λ denotes the position of the source. α is a unit 3D vector (∈ S 2 ) that gives the direction of one X-ray. Note here that the system is perfectly aligned: the v-axis is parallel to the rotation y-axis (θ = η = 0). The u-axis is perpendicular to the direction of the source (φ = 0) and the detector is not shifted (u 0 = v 0 = 0).

independent coefficients of the 3 × 4 projection matrix in homogeneous coordinates (see e.g. [START_REF] Hartley | Multiple View Geometry in Computer Vision Second Edition[END_REF]). If the detector rows and columns are known to be perpendicular with the same sampling in both directions (i.e. square pixels) then two degrees of freedom are eliminated and the task reduces to estimating nine projection-specific geometric parameters (three for the source position, three for the detector position and another three for the orientation of the gantry). As described in the previous section, the circular trajectory we are considering in this work is completely described by 7 global geometric parameters. The only projection-specific parameter is the rotation angle.

Imaging-based calibration methods fall into two broad categories. One category consists of the off-line methods [START_REF] Cho | Accurate technique for complete geometric calibration of cone-beam computed tomography systems[END_REF][START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF][START_REF] Gullberg | Estimation of geometrical parameters and collimator evaluation for cone beam tomography[END_REF][START_REF] Mennessier | Direct determination of geometric alignment parameters for cone-beam scanners[END_REF][START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF][START_REF] Rizo | Geometric calibration method for multiple-head cone-beam SPECT system[END_REF][START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF]: they all require pre-scanning of a calibration phantom, usually made up of small ball bearings (BB) whose relative positions are accurately known. Then the theoretical projections of the BBs (which depend on the geometric parameters) are compared with their actual projections to derive -iteratively [START_REF] Gullberg | Estimation of geometrical parameters and collimator evaluation for cone beam tomography[END_REF][START_REF] Rizo | Geometric calibration method for multiple-head cone-beam SPECT system[END_REF] or analytically [START_REF] Cho | Accurate technique for complete geometric calibration of cone-beam computed tomography systems[END_REF][START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF][START_REF] Mennessier | Direct determination of geometric alignment parameters for cone-beam scanners[END_REF][START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF][START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF] -the calibration parameters. In [START_REF] Mennessier | Direct determination of geometric alignment parameters for cone-beam scanners[END_REF], they solve the complete calibration problem and analytically derive all nine parameters for each projection.

The other group of techniques consists of on-line techniques. All methods in this category solve the calibration problem without a specific cali-bration scan of a calibration object. They only use the data from the projections of the imaged object. Beyond this common feature, this group encompasses substantially different techniques. In [START_REF] Kingston | Reliable automatic alignment of tomographic projection data by passive auto-focus[END_REF][START_REF] Kyriakou | Simultaneous misalignment correction for approximate circular conebeam computed tomography[END_REF][START_REF] Muders | Stable and Robust Geometric Self-Calibration for Cone-Beam CT Using Mutual Information[END_REF], they minimize a cost function, whose evaluation requires the reconstruction of the object from the current estimate of the geometric parameters. The metric is based on entropy in [START_REF] Kyriakou | Simultaneous misalignment correction for approximate circular conebeam computed tomography[END_REF], the L 2 -norm of the image gradient in [START_REF] Kingston | Reliable automatic alignment of tomographic projection data by passive auto-focus[END_REF] or the mutual information between re-projected image and projection data in [START_REF] Muders | Stable and Robust Geometric Self-Calibration for Cone-Beam CT Using Mutual Information[END_REF]. The limitation of such methods is the computational load, which may not fit clinical workflow (though [START_REF] Kingston | Reliable automatic alignment of tomographic projection data by passive auto-focus[END_REF] limits this drawback by only reconstructing a fraction of the volume). In [START_REF] Ouadah | Self-calibration of cone-beam CT geometry using 3D-2D image registration[END_REF], they use the 3D reconstruction of a planning CT and compute projection-specific geometric parameters by registering the actual projections with the re-projected CT image. Other works in this category utilize the redundancy of the projection data (ie. the DCCs). In the 2D parallel beam case, Basu and Bresler [START_REF] Basu | Feasibility of tomography with unknown view angles[END_REF][START_REF] Basu | Uniqueness of tomography with unknown view angles[END_REF] solve uniquely and efficiently the problem of unknown projection angles and shifts with the Helgason-Ludwig DCCs. Some works use the trivial "opposite-ray" condition [START_REF] Meng | Online Geometric Calibration of Cone-Beam Computed Tomography for Arbitrary Imaging Objects[END_REF][START_REF] Panetta | An optimizationbased method for geometrical calibration in cone-beam CT without dedicated phantoms[END_REF][START_REF] Patel | Self-calibration of a cone-beam micro-CT system[END_REF]. In [START_REF] Panetta | An optimizationbased method for geometrical calibration in cone-beam CT without dedicated phantoms[END_REF], this DCC, which normally only applies in the central plane (the plane of the trajectory), is extended to cone-beam projections of a particular class of symmetric 3D object functions and shows accurate calibration results when approximated in a central region of a generic object. More closely related to our work is a series of publications on epipolar consistency conditions [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF][START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF][START_REF] Maass | Geometrical Jitter correction in computed tomography[END_REF]. These DCCs are based on the Grangeat theorem and relate the derivative of the 2D Radon transform of the projections to the derivative of the 3D Radon transform of the imaged object.

The comparison of previous works is not easy due to the parametrization which may differ with authors. Nevertheless, it is widely documented that the detector shift u 0 and the in-plane angle η are of crucial importance [START_REF] Daly | Geometric calibration of a mobile C-arm for intraoperative cone-beam CT[END_REF][START_REF] Noo | Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography[END_REF][START_REF] Smekal | Geometric misalignment and calibration in cone-beam tomography[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF]. On the other hand, [START_REF] Muders | Stable and Robust Geometric Self-Calibration for Cone-Beam CT Using Mutual Information[END_REF][START_REF] Yang | A geometric calibration method for cone beam CT systems[END_REF] demonstrated that the two outof-plane angles (φ and θ) may be set to zero without affecting the image quality if their true values are kept below 2 • (which is a reasonable manufacturing accuracy requirement). Finally, miscalibration of the source-tocenter and source-to-detector distances does not introduce artifacts in the reconstructed volume and are therefore not calibrated. However, these two parameters affect the magnification of the reconstructed volume, which would not be acceptable in some cases, such as a metrology-oriented application.

M E T H O D S

Cone-beam DCCs for a linear trajectory

Data consistency conditions are conditions which must be satisfied by the projection data in order for them to be the image of an object function through the forward projection model described in Equation 188. DCCs have been applied to various CT artefact correction techniques, e.g. mo-tion compensation [START_REF] Clackdoyle | Fanbeam data Consistency Conditions for applications to motion detection[END_REF][START_REF] Frysch | Rigid motion compensation in interventional C-arm CT using consistency measure on projection data[END_REF][START_REF] Yu | Data consistency based translational motion reduction in fan-beam CT[END_REF] and beam hardening correction [START_REF] Tang | Data consistency condition-based beam-hardening correction[END_REF]. The simplest condition is the "opposite-ray condition". In parallel projection geometry, it states that the projections must be even: R f ( α, s) = R f (-α, -s), where R f denotes the 2D Radon transform of an object function f . This condition was applied to the calibration problem in [START_REF] Panetta | An optimizationbased method for geometrical calibration in cone-beam CT without dedicated phantoms[END_REF]. Still in the parallel geometry, there exists a complete set of DCCs, known as Helgason-Ludwig DCCs (see [START_REF] Helgason | The Radon transform. Progress in mathematics[END_REF][START_REF] Ludwig | The radon transform on euclidean space[END_REF] and standard textbooks on the Radon transform, e.g. [START_REF] Deans | The Radon Transform and some of its applications[END_REF][START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF]), which relates the n-th order moments of each projection to a homogeneous polynomial of order n. In the 2D divergent geometry (fanbeam projections), complete DCCs, similar to Helgason-Ludwig polynomial conditions, were derived in [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF] for the particular case of an X-ray source moving along a line. We will be using the order-0 case, which was known much earlier than the latter work (see their various guises in [START_REF] Chen | A new data consistency condition for fanbeam projection data[END_REF][START_REF] Finch | Sums of homogeneous functions and the range of the divergent beam x-ray transform[END_REF][START_REF] Noo | Image reconstruction from fan-beam projections on less than a short scan[END_REF][START_REF] Wei | Integral invariants for computed tomography[END_REF] and [START_REF] Tang | The mathematical equivalence of consistency conditions in the divergent-beam computed tomography[END_REF] for a review).

The description of cone-beam pair-wise consistency conditions follows [START_REF] Levine | Consistency Conditions for Cone-Beam CT Data Acquired with a Straight-Line Source Trajectory[END_REF]. Let s λ i and s λ j be two source positions and L i , j be the line connecting them. Suppose that both projections are acquired with one common flat detector, parallel to L i , j . Any plane containing L i , j intersects -if it doesthe detector on a row, parallel to L i , j , which we will index with k. We will denote that plane P i , j ,k . The situation in P i , j ,k reduces to a pair of fanbeam projections along the virtual linear trajectory L i , j and with the k-th detector row playing the role of the 1D fanbeam detector. The order-0 DCCs state that:

Lemma 4. For any pair of projection indices i , j and any integer k, let:

G i , j ,k = π 2 -π 2 g λ i , α k φ cos φ d φ, ( 189
)
where α k φ denotes a unit vector in P i , j ,k , φ denotes the angle between α k φ and the perpendicular to L i , j in that plane. Furthermore, the line L i , j is assumed to not intersect the support of the object function f . See Fig 32 and 33. If the data are consistent, then:

G i , j ,k -G j ,i ,k = 0, (190)
Let c i , j ,k denote the square difference of the LHS of Equation 190. The sum C i , j = k c i , j ,k is a measure of the pair-wise consistency between two cone-beam projections g (λ i , •) and g (λ j , •).

Re-sampling in a virtual detector

These DCCs only apply if the detector is parallel to the virtual linear trajectory L i , j connecting two source positions. In the circular trajectory we are considering in this work, this detector condition is obviously not fulfilled. To remedy this problem, each pair of projections is re-sampled onto a virtual detector V i , j by means of a backprojection. The virtual detector is placed at the origin of the world system of coordinates and oriented in such a way that the rows and columns of V i , j are parallel to L i , j and the axis of rotation, respectively. The situation is illustrated in Figure 33.

The orientation of the virtual detector allows a simple evaluation of the integral G i , j ,k in Equation (189) by changing the φ-variable to the u-pixel coordinate of the virtual detector with:

u = v 2 k + D 2 virt tan φ,
where v k is the intercept of the plane P i , j ,k with the virtual detector's vaxis and D virt denotes the distance from the source to the virtual detector. Applying this change of variables leads to:

G i , j ,k = 1 v 2 k + D 2 virt R g (λ i , u) v 2 k + D 2 virt u 2 + v 2 k + D 2 d u. ( 191
)
Note that the weight inside the integral is exactly cos φ. The change of variables has moved this cosine term from the denominator to the numerator.

When applying these DCCs to the calibration problem, we note that the backprojection onto the virtual detector will use the projection geometry as input. Hence the dependency of the cost function (described in the next section) on the calibration parameters via this backprojection.

Figure 33: Two sources on a circular trajectory. Both projections are backprojected in a virtual detector, parallel to the line connecting s λ i and s λ j .

The consistency metric

Estimation of the geometric calibration parameters is achieved by minimizing a cost function based on the pairwise consistency conditions described above. Let p = (φ, θ, η, u 0 , v 0 , R, D) denote the 7-uple of sought parameters (remember that the rotation angle λ is accurately known). We define the cost function C (p) as follows:

C (p) = (i , j )∈Ω C i , j = (i , j )∈Ω k c i , j ,k , (192)
where C i , j was defined above with the dependence on p buried in the backprojection onto the virtual detector, Ω is the chosen subset of pairs of projections to which the DCCs are applied. The size of the virtual detector V i , j is computed to account for the distortion resulting from the backprojection step (see details is Section 2.6.3.1). For each pair of projections, the cost function is evaluated over all rows k of V i , j .

The computation of the cost function can be summarized in the following algorithm: Backproject projections onto V i , j .

5:

Pre-weight the virtual projections acc. to Equation 191. end for 12: end procedure

N U M E R I C A L E X P E R I M E N T S O N S I M U L AT E D D ATA

We first studied the properties of our cost function on simulated projections of a Shepp-Logan phantom and estimated the accuracy that can be expected from our method. All simulated projection data were generated with the RTK software package [START_REF] Rit | The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK)[END_REF]. All reconstructions were computed with the FDK algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] available in RTK.

In all our experiments, the set Ω was composed of 27 projections pairs, constructed as follows: nine equally spaced projections (spaced by 40 • ) were selected and all possible pairs were included in Ω, except those separated by ±160 • . This particular choice for Ω arose from a trade-off between the computational load and the amount of data we inject in the cost function for robust parameter estimation. Pairs separated by ±160 • were removed because they are too close to the limit situation where the line L i , j (hence the virtual detector) would be perpendicular to the physical detectors. Also, the maximum separation of the remaining pairs was 120 • which eliminated any risk of the connecting line intersecting the scanned object.

Cost function study

We first studied the behaviour of our cost function on the simulated projections of a 3D Shepp-Logan phantom [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. Projections were simulated with a perfectly aligned system (R = 100, D = 200 and all other geometric parameters set to zero). Then, we computed the cost function as a 1D-function of each separate parameter, over a symmetric interval [-2, 2] (in degrees for η, θ and φ and in millimeters for u 0 , v 0 , R and D). Figure 34 shows corresponding plots. Note that the cost function has very low dependence on the two distances R and D. For this reason, these two parameters will not be optimized in our calibration method. Our procedure focuses on the remaining 5 global parameters. The plots in Figure 34 indicate that the cost function is locally convex with respect to each of them. The convexity of the 1D-functions does not guarantee the convexity of the multi-dimensional cost function but is still encouraging for the optimization procedure to find a suitable minimum. Of the five parameters, the vertical shift v 0 shows the least sensitivity to the DCCs. This fact has a direct effect on the errors we obtained with simulated data (see section below).

Calibration on simulated data

We applied our calibration method to simulated projections of a 3D Shepp-Logan phantom. The data were generated using the misalignment parameter values indicated in the first row of Table 3. The simulated projections were 256 × 256 pixels, with pixel size set to 0.25mm. The source-to-center and source-to-detector distances were set to 100 mm and 160 mm respec- tively. The larger half-length of the outer ellipsoïd of the Shepp-Logan phantom was 15 mm. We used the Numpy implementation of the order-0 minimization method from Powell [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF]. Results are recorded in Table 3.

In the first experiment (see Exp. 1 in Table 3), the calibration procedure was initialized with a nominal geometry (all five parameters equal to 0). This initialization corresponds to the best guess we could make on the real µCT system, which is designed to be perfectly aligned. In Figure 35, we present one frontal slice of the 3D numerical Shepp Logan phantom (top left). We reconstructed the 3D volume from the simulated projections data using three different geometries: 1) the true geometry, 2) the DCCcalibrated geometry (resulting from our minimization procedure) and 3) the nominal geometry. The reconstructions are shown in Figure 35. The reconstruction with the nominal geometry shows severe artifacts (see middle left in Figure 35) with a RMSE of 0.311 when compared to the 3D numerical phantom (top left). No visible difference between the two reconstructions with the true geometry (top right) and with DCC-calibrated geometry (middle right) is apparent. In both cases, the quality of the reconstruction is significantly improved, with RMSE of 0.107 and 0.106 respectively. Note also that the procedure can easily be extended to a short scan trajectory. In the bottom row of Figure 35, we present a 220 • short-scan reconstructions with the nominal geometry (left) and the DCC-calibrated geometry The second experiment focuses on the dependency of the cost function on the initial guess. We ran the procedure with 100 random initial values taken from a normal distribution and computed the mean and standard deviation of each geometric parameters. The results are presented in Table 3, Exp. 2. We noticed that the out-of-plane angle θ and the vertical shift of the detector v 0 deviated significantly from their true values (mean errors: 1.28 and 2.16 respectively), with large variability (standard deviations: 1.48 and 4.11 respectively). But this does not affect the quality of the reconstruction as shown in Figure 35. We observed that the imaged object was well reconstructed but possibly at a slightly different place in space and the reconstruction with DCC-calibrated geometry had to be manually registered on the reference Shepp-Logan by a small translation in the v direction. On the other hand, when one of these two parameters was fixed to its true value, the optimized value of the other was estimated with a high precision and with very low variance (see Table 3, Exp. 3 and Exp. 4). This reveals the difficulty of evaluating the quality of a calibration procedure. To better illustrate this compensation phenomenon, we computed the cost function value as a function of θ and v 0 (the other 3 parameters being fixed to their true values). The plot in Figure 36 reveals a long flat valley in a direction which is a linear (for small θ) combination of θ and v 0 directions. All the values of θ and v 0 along this valley minimize the cost function and visual inspection of the resulting reconstruction suggested equally good quality.

N U M E R I

C A L E X P E R I M E N T S O N O U R µ -C T S Y S T E M

Description of the experimental set-up

We applied our method to real data acquired on the micro-CT system depicted in Figure 29. The X-ray source is fixed. The imaging object is placed on a turntable. The distance between the source and the turntable can be manually adjusted. The source to detector distance can also be manually adjusted, but both distances remain fixed during one complete 360 • scan. The detector size is 2048 × 2048 pixels with pixel size of 17.09 2 µm 2 . The X-ray tube voltage was set to 30 kV, the current was set to 0.25 mA and exposure time was set to 2 s. To reduce the computational burden, the projections were down-sampled by a factor of 8 to 256 × 256 pixels.

In each experiment, 360 equally-spaced projections were acquired. Darkfield and flat-field corrections were applied to raw-data. The negative-log transform was then applied so that the projection data correspond to the line-integral model described in Equation 188. We report results on three different datasets: the first one is the projection data of the plastic cap from a tube of glue, approximately 1 centimeter wide. The rotating support plat-form was in the flat-field images and therefore subtracted from the projection of the glue-cap. The glue-cap was small enough to be completely contained in the projections. Consequently, no projections were truncated. The values of R and D were physically measured to be 219 mm and 295 mm. The second dataset was acquired from a sample of concrete foam. All projections were truncated in the direction of the rotation axis (axial truncation). Measures of R and D were 114 mm and 137 mm respectively. For this sample, a 0.4 mm aluminum filter was placed in front of the X-ray tube to harden the X-ray beam to make the projection data better fit the line integral model. The third data consisted in a piece of sponge placed into a plastic syringe. Projections were also axially truncated. Measures of R and D were 195 mm and 259 mm respectively. See pictures and sample projections of the three objects in Figure 37.

Reconstruction with complete data

The calibration method was applied to the projections acquired from each scan, using the nominal geometry as first guess. The output values are indicated in Table 4. Each of the eight row in this table corresponds to a different subset of projections, from which the cost function was computed.

The first one was composed of nine equally spaced projections, starting with projections at angle 0. Each subsequent subset was shifted by 5 projections (5 degrees). Figure 38 shows coronal and transverse slices from the reconstructed images with the nominal geometry and compares them to reconstructions with a geometry estimated using an off-line marker-based method and our DCC-calibrated geometry. The alignment problem described in Section 2.5.2 was encountered here too and the two calibrated reconstructions were registered manually in the y direction for comparison. Note first that sub-degree angular misalignments and sub-millimeter detector shifts lead to severe artifacts in the reconstruction, especially at 38). Second, the image quality was significantly improved when reconstruction was computed with the DCC-calibrated geometry. The edges are sharp as illustrated by the profiles in Figure 38. Of course, the calibration procedure does not correct for other CT artifacts which degrade both un-calibrated and calibrated reconstructed images (e.g. cupping, probably due to beam hardening, and ring artifacts).

Reconstruction with axially truncated data

This section explains how our calibration procedure can deal with axially truncated data with application to the truncated data acquired on the same µ-CT system (Figure 37 middle and right).

Handling axial truncation

Our cost function is the sum of square differences between integral over rows of the virtual detector. For that reason, truncation in the v-direction does not cause any difficulty as long as there is no truncation in the udirection. This feature is specific to the nature of the DCCs used in the cost function. In our implementation, care must be taken at the backprojection level because the square physical detector is backprojected to a trapezoidal shape on the virtual detector, with horizontal pixel rows backprojected to oblique pixel rows of varying angle (except for the central line, which remains horizontal). The situation is depicted in Figure 39. The virtual projection can therefore be limited to those horizontal rows of the virtual detector that are not truncated (rows between the two dashed lines on Figure 39 right).

Results

The calibration procedure was applied to the concrete and the sponge datasets. The nominal geometry served as initial guess. For the concrete sample, the scanning distances R and D were set to 114 mm and 137 mm respectively. The resulting cone-angle was approximately 14 • . For the sponge sample, R = 195 mm and D = 259 mm. Axial and transverse slices of the reconstructed volumes are shown in Figure 40 and 41. In the uncalibrated reconstructions, small structures of the object are barely distinguishable. In the calibrated reconstruction of the concrete sample, though cone-beam and beam-hardening artifacts are still present, the detailed structures (air bubbles in the concrete foam) are much more sharply reconstructed.

C O N C L U S I O

N A N D F U T U R E W O R K
We proposed an on-line calibration method to estimate five geometric parameters of a µ-CT system. The method is based only on consistency of the "production" scan. It requires no prior (off-line) calibration scan. The quality of the reconstructed images in the experiments compares with the robust "classical" marker-based calibration method. Furthermore, the calibration method can correctly handle axially-truncated data, which is an untypical feature for DCC-based application. The design of our cost function can probably be refined. A short study on the individual contribution of a pair of projections revealed that pairs angularly separated by more than 90 degrees contributed more than close pairs. Hence, a cost function built from such pairs may convey more independent information and hence lead to more robust estimation. Another question is related to the dependency of the cost function on the object. We have carried out some simulations (similar to those in Figure 34 and Figure 36) on objects with sharp edges (a simplex-like simulated object) or plate-like objects (very small extent in the v-direction). In all cases, the cost function behaved similarly to the Shepp-Logan study, with regards to each individual parameter or with regards to the (θ, v 0 ) pair. However, the cost function behaved differently when the plate-like object was placed in the central plane (containing the source trajectory). But, in this case, the geometry collapses to fanbeam, with its own geometric parameters (for example, θ plays no role).

The investigation of the interplay between geometric parameters is a possible future direction of research. Figure 36 reveals that a large error on one parameter can be compensated by a large error on the second in terms of consistency. We are also extending this work in two directions. The first one applies the same principles to estimate projection-specific calibration parameters, by using a similar cost function for each projection. Second, the comparison of our method with the work in [START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF], later described as Epipolar Consistency Conditions [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF]. ECCs are also applied to pairs of projections and use a similar geometry of lines on the two detectors (as shown in Figure 33). However, the theoretical foundations are different because the ECCs are based on Grangeat's formula and require the computation of a derivative. Whether these conditions are equivalent to the conditions used in our work still needs to be understood and is ongoing work.

I N V E S T I G AT I O N O F T H E V A L L E Y -S H A P E D C O S T F U N C T I O N

This section is an addendum. It was not published in the original paper [START_REF] Lesaint | Calibration for Circular Cone-Beam CT Based on Consistency Conditions[END_REF]. We investigate the behaviour of the cost function in the 2D plane of the parameters θ and v 0 . Figure 36 revealed a flat valley around the nominal position (θ, v 0 ) = (0.2, 0.5). At any point along that valley, the consistency is minimal. In other words, an error on -say -θ is compensated (in terms of consistency) by an error on v 0 . More precisely, we prove that the consistency cost function is insensitive to small (to be clarified) errors on θ.

To draw Figure 36, we sampled the cost function C over the parameters θ and v 0 for θ ∈ 0.2 + [-2, 2] and v 0 = 0.5 + [-2, 2]. For each v i 0 , we denoted

θ i = argmin θ C (θ, v i 0 ). ( 193
)
This value was just extracted from the samples of C (no minimization procedure was launched). The obtained points lie along the valley. We then fit these points by linear regression. The equation we obtained was

v 0 ≈ 159.352 × θ -0.064, (194)
with very strong correlation (r > 0.999). This linear relation must be true only for small values of θ since we anticipate a trigonometric relation between these two variables. A study of the cost function along the valley for larger values of θ (θ ∈ 0.5±10 • ) also shows that this compensation does not occur for too large errors (see Figure 42).

To derive the geometric relation between θ and v 0 , we consider (without loss of generality) one single projection at rotation angle λ = 0, with a detector aligned in the x-direction, i. e. η = φ = 0 = u 0 . The focal distance is f . The nominal values of θ and v 0 are 0.2 • and 0.5 mm respectively. When changing the value of θ without affecting the other parameters (θ denotes the modified value), the detector is rotated around the position of the source (see Figure 43). This rotation introduces inconsistency. E. g. , the ray from the source to the detector's centre is significantly displaced. We define v 0 to be the detector's offset such that the ray through the displaced detector's centre coincide with the same ray in the nominal position. We have

v 0 = f tan θ, ( 195
) v 0 = f tan θ , (196)
so that

v 0 = f tan θ + v 0 -f tan θ. ( 197
)
With the numerical values of our simulations ( f = 160mm, v 0 = 0.5 and

θ = 0.2 • ) one gets v 0 = 160 tan θ -0.0585, (198)
which is close to Equation 194 (with the approximation tan θ ≈ θ valid for small values of θ ). This short study illustrates that the consistency-based cost function is insensitive to small variations of the out-of-plane angle θ (at least for errors of less than 0.7 • in amplitude, see Figure 36), provided that the detector is offset accordingly to realign at least one ray. If the out-of-plane angle is modified to θ , the ray from s λ to the centre of the detector is displaced accordingly, introducing inconsistency. The inconsistency is counteracted by adjusting v 0 so that the ray returns to its original position.

II3 F B C C A N D G C C : A C O M PA R AT I V E S T U D Y
In this chapter, we carry out a comparative study between FBCC presented in Chapter II-2 and the Grangeat Consistency Conditions (GCC) first introduced in [START_REF] Debbeler | A New CT Rawdata Redundancy Measure applied to Automated Misalignment Correction[END_REF] 30 . The main result of this chapter establishes the intrinsic 30 For a reason that will become clear later in this chapter, the name of Epipolar Consistency Conditions, used in [START_REF] Aichert | Epipolar Consistency in Transmission Imaging[END_REF] is not used here. We prefer to use Grangeat Consistency Conditions to emphasize that the DCC essentially rely on Grangeat's result.

In fact, FBCC are as epipolar as the GCC since the word epipolar only refers to the geometry of lines in the detectors.

In both FBCC and GCC, the same set of epipolar lines are incorporated in the conditions.

connection between the GCC and FBCC (Proposition 16). By intrinsic, we mean that this relation is independent of any detector. Then, for practical implementation of both conditions, we establish the same result for pairs of source and detector in general positions. To this end, we derive the implementation of FBCC in the physical detector (Proposition 18) and establish an epipolar version of Grangeat's result (Proposition 20). We then verify that the intrinsic relation is also satisfied by the expressions in the physical detector. Finally, we illustrate this relationship with numerical simulations. The reading of this chapter may be facilitated by the prior reading of the Appendix 1.3, which introduces the concepts of the projective geometry, homogeneous coordinates and homography.

T H E E P I P

O L A R G E O M E T R Y O F O N E P A I R O F P R O J E C T I O N S
We consider an object function µ whose support contains the origin O of the PCS. The source trajectory never intersects the support of µ (in particular, the source is never in O). Let g λ 1 and g λ 2 be a pair of projections. The source positions are denoted s 1 and s 2 , with coordinates in the PCS (x λ 1 , y λ 1 , z λ 1 ) and (x λ 2 , y λ 2 , z λ 2 ) respectively. The line connecting these two points is called the baseline and is denoted b. Two cone-beam projections define an epipolar coordinate system (ECS) (O, e e

x , e e y , e e z ) illustrated in Figure 44 and defined by

• e e x = s 2 -s 1 s 2 -s 1
is the direction of the baseline b,

• e e y = s 1 × s 2 s 1 × s 2 is the unit normal to the plane containing both source positions and the origin O,

• e e z = e e x × e e y is such that the ECS is a right-handed coordinate system.

We let R 12 = e e x e e y e e z denote the rotation matrix from ECS coordinates to PCS coordinates and d > 0 denotes the distance from the origin O to the baseline b. Note that in the ECS, every point on the baseline has coordinates (x, 0, d ) for some x ∈ R. In particular, we define the scalars α 1 and α 2 so that the source positions s 1 and s 2 have coordinates (α 1 , 0, d ) and (α 2 , 0, d ) respectively in the ECS (Figure 45). Note that we have We also define spherical coordinates (θ, φ) so that a unit vector α ∈ S 2 has coordinates, in the ECS (see Figure 45), α(θ, φ) = (sin φ, sin θ cos φ, -cos θ cos φ).

(x λ i , y λ i , z λ i ) = R 12 (α i , 0, d ).
(199)

We call epipolar plane any plane which contains the baseline b. The 1D family of epipolar planes can be parametrized by the angle θ ∈ ]-π/2, π/2] that a particular plane makes with the plane (O, e e

x , e e z ). We let E θ denote an epipolar plane. Each of these planes has a normal unit vector denoted β(θ), whose coordinates in the ECS are β(θ) = (0, cos θ, sin θ).

(

A plane E θ intersects each physical detector in one epipolar line, denoted l θ,i (i = 1, 2). All epipolar lines l θ,i in one detector meet in one point, called the epipole and denoted ε i . The epipole ε 1 is the projection of s 2 through the source s 1 (and conversely for ε 2 ). The situation is depicted in Figure 46.

3.2 A D E T E C T O R -L E S S V E R S I O N O F P A I R -W I S E D C C

Grangeat Consistency Conditions (GCC)

In this section, we briefly recall Grangeat Consistency Conditions. The theorem of Grangeat was introduced in Section 1.7 with the use of the inter- mediate function G g , defined for any projection g λ and for all β ∈ S 2 by

s 1 s 2 ε 1 ε 2 E θ l θ,1 l θ,2
G g λ, β = S 2 δ α • β g λ α d α. ( 201
)
Grangeat theorem relates G g λ, β to the 3D Radon transform of the object function µ

G g λ, β = ∂ ∂s Rµ β, s s λ • β . ( 202
)
Let us consider an epipolar plane E θ which contains both sources s 1 and

s 2 . Since s λ 1 • β(θ) = s λ 2 • β(θ), one has the Grangeat consistency condition S 2 δ α • β(θ) g λ 1 α d α = S 2 δ α • β(θ) g λ 2 α d α. ( 203
)
Of course, in some particular cases, a plane E θ may contain more than two source positions. This is the case in Chapter II-1: the source trajectory is assumed planar. The plane E 0 (it is denoted H e y ,0 in Chapter II-1) contains all source positions, so the equality of Equation 203 propagates to the whole trajectory. It is also the case when the trajectory of the source is linear as in [START_REF] Clackdoyle | Necessary and Sufficient Consistency Conditions for Fanbeam Projections Along a Line[END_REF] and in [START_REF] Lesaint | GCC and FBCC for linear tomosynthesis[END_REF], in the case of linear tomosynthesis.

Fan-Beam Consistency Conditions (FBCC)

The second set of DCC that we are considering in this chapter is the FBCC, which were applied to the calibration problem in Chapter II-2. With the spherical coordinates defined in Equation 199 and for a particular plane E θ , the pairwise FBCC can be formulated as

π 2 -π 2 g λ 1 α(θ, φ) cos φ dφ = π 2 -π 2 g λ 2 α(θ, φ) cos φ dφ. ( 204
)
We denote G f λ, β(θ) the integral involved in this equation.

The intrinsic relation between G g and G f

We can now establish the main result of this section, which relates FBCC and GCC.

Proposition 16. For all θ ∈ -π 2 , π 2 , G g λ, β( θ) = ∂ ∂θ G f λ, β(θ) θ= θ. (205)
The relation in Equation 205 is independent of the detectors. Saying that the pair of projections satisfy the FBCC means that the functions G f (λ 1 , •) and G f (λ 2 , •), as functions of the angle θ, are equal for all θ. Their derivatives (with respect to θ) must also be equal and the projections satisfy the GCC. We proved the

Corollary 17. If a pair of projections satisfies the FBCC, then it must also satisfy the GCC.

To the best of our knowledge, this result is new.

Proof of Proposition 16. We start with the expression of G g in terms of the spherical coordinates (θ, φ). We obtain

G g λ, β( θ) = + π 2 -π 2 +π -π δ cos φ sin(θ -θ) g λ α(θ, φ) cos φ dθ dφ. (206)
Again, the homogeneity property of δ plays an important role by moving the cos φ from the numerator to the denominator, yielding

G g λ, β( θ) = +π -π δ (sin(θ -θ)) + π 2 -π 2 g λ ( α(θ, φ)) cos φ dφ dθ (207) = ∂ ∂θ G f λ, β(θ) θ= θ, (208)
where we obtained the last line by applying the Proposition 24 of Appendix 1.1 (page 124).

The result of Proposition 16 is of theoretical interest. It states that the FBCC are stronger conditions than GCC, in the sense that if FBCC are satisfied by a pair of projections, the GCC are necessarily also satisfied. Whether they are strictly stronger (i. e. we could find a pair of projections which are consistent with respect to the GCC but inconsistent with respect to FBCC) is still unclear. It is true in the particular tomosynthesis geometry. We derived such an example in [START_REF] Lesaint | GCC and FBCC for linear tomosynthesis[END_REF]. In view of applications, we will now establish the same relation in terms of the detector coordinates. To do so, we need first to derive a formula for G f on the detector. This is the purpose of the next section.

F A N -B E A M D C C I N T H E P H Y S I C A L D E T E C T O R

In this section, we derive the algorithm to compute the FBCC in the physical detector. The situation is as above except that we now add physical detectors, whose positions and orientations are entirely encoded in the projection matrices P 1 and P 2 . Both projections are assumed in very general position. We only assume that the baseline b does not intersect the support of the object 31 . In particular, the baseline never passes through the origin 31 This assumption is not too restrictive with respect to a circular trajectory for instance. If projections are acquired at a reasonable sampling rate (say every five degrees at least), there will be sufficiently many pairs which fulfil the condition.

O. We also assume that neither detector is parallel to the baseline.

From now until the end of this chapter, we focus on one projection only and drop the subscript i = 1 or 2 for ease of notation. Of course, the epipolar geometry of the pair is still ubiquitous, through the rotation matrix R 12 and the ECS. For a particular epipolar plane E θ making an angle θ with the reference plane (O,e e x , e e z ) 32 , the FBCC require the computation of the 32 Note that θ is also the angle between β(θ) and e e y .

function G f defined by

G f (θ) = G f λ, β(θ) = π 2 -π 2 g λ α(θ, φ) cos φ dφ. (209)
In Chapter II-2, the computation of G f is implemented on a virtual detector whose axes coincide with the e e

x and e e y directions of the ECS. The epipolar plane E θ intersects the virtual detector in a line parallel to the uaxis with a vertical offset which we denote v θ . Note that we have the re-lation v θ = d tan θ. By changing the φ variable in Equation 209 for the u variable (see Figure 45), with

u = α + d cos θ tan φ (210) dφ = v 2 θ + d 2 (u -α) 2 + v 2 θ + d 2 du (211) cos φ = v 2 θ + d 2 (u -α) 2 + v 2 θ + d 2 , ( 212
)
we obtain

G f (θ) = +∞ -∞ g λ (u, v θ ) (u -α) 2 + v 2 θ + d 2 du. ( 213
)
The integral involved in Equation 213 is a 1D integral over a line parallel to the u-axis of the virtual detector. When projected onto the physical detector, the line of integration will not be parallel to one of the axes (see Figure 47). To address this issue, we turn the 1D integral into a 2D integral, using the Dirac δ distribution

G f (θ) = δ(v -v θ ) g λ (u, v) (u -α) 2 + v 2 + d 2 du dv. (214)
In view of implementing this function in the physical detector, we will change the (u, v) variables of the virtual detector to the coordinates of the physical detector, which we denote with capital letters (U ,V ). The final result of this change of variables is Proposition 18 (FBCC in physical detector).

G f (θ) = 1 cos e o z , e e x l θ f (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 g λ (U ,V ) (U -U ε ) 2 + (V -V ε ) 2 dσ ( 215
)
where l θ is the epipolar line and (U ε ,V ε ) are the coordinates of the epipole ε.

The integral over l θ must be understood in the sense of curve integral 33 . 33 Let C be a differentiable curve parametrized with

x(t ) for t ∈ T . We have C f ( x) dσ = T f ( x(t )) | x (t )| dt . R E M A R K S :
The term before the integral is the cosine of the angle made by the normal of the physical detector and the baseline. It does not depend on (U ,V ) (but it is projection dependent, because e o z is). The first weighting term inside the integral is the usual cosine of the incidence angle. Note that in the virtual detector, the projection values were also weighted by the cosine of the incidence angle (up to the multiplicative constant d ), but this incidence angle was with respect to the virtual detector. Shifting to the physical detector substitutes the cosine of the virtual incidence angle to the cosine of the physical incidence angle. The second weighting term is

U ε i 2 + V -V ε i 2 (216)
is actually the distance from the pixel to the baseline b. It corresponds to Equation [START_REF] Cho | Accurate technique for complete geometric calibration of cone-beam computed tomography systems[END_REF] in [START_REF] Clackdoyle | Fanbeam data Consistency Conditions for applications to motion detection[END_REF]. The correspondence is not immediate but the denominator in their equation is precisely the distance of the pixel to the line connecting the source positions.

Proof of Proposition 18. The rest of this section is dedicated to the proof of the Proposition 18. It consists in the following four steps:

• Compute the projection matrices from ECS coordinates to the virtual (u, v) coordinates and physical (U ,V ) coordinates.

• Derive the homography that leads to the change of variables.

• Change the variables.

• Properly handle the composition of the Dirac δ distribution with the change of variables.

At every stage of the proof, we try to work out the geometric meaning of the expressions, instead of simply deriving weighting terms which would obscurely depend on the parameters of an homography H (see below). Most calculations will be done in the ECS. The first step consists in deriving the projection matrix from ECS coordinates to the virtual and physical DCS.

S T E P 1 :

P R O J E C T I O N M AT R I C E S F R O M E C S C O O R D I N AT E S T O V I R - T U A L A N D P H Y S I C A L D E T E C T O R S .
The position and orientation of the source and the (physical) detector is encoded in the projection matrix P , which can be be written in several forms (see Chapter I-2):

P = K R I 3 |-s = A| -A s with K =     -f 0 u 0 0 -f v 0 0 0 1     (217) R =     -e o x - -e o y - -e o z -     (218) s =     x λ y λ z λ     (219) A = K R. (220)
The matrix P converts the 3D coordinates in the PCS to the 2D coordinates in the (physical) detector. We let P p denote the projection matrix from the ECS coordinates to the physical detector coordinates. We let R12 be the rotation matrix from the ECS to the PCS, in homogeneous coordinates.

R12 =        | | | 0 e e
x e e y e e z 0

| | | 0 0 0 0 1        = R 12 0 0 1 . (221)
The matrix P p is obtained by right-multiplying P by R12 .

P p = P R12 (222) = K R I 3 | -s R 12 0 0 1 (223) = K R R 12 | -s (224) = K RR 12 I 3 | -R T 12 s . (225)
The matrix P p can be decomposed as

P p = K p R p I 3 |-s p = A p |-A p s p with K p = K (226) R p = RR 12 (227) A p = K RR 12 ( 228
)
s p = R T 12 s. ( 229
)
We now derive the projection matrices from ECS coordinates to the virtual detector, denoted P v . We built the matrix from the geometrical description of the source and virtual detector (see Chapter I-2). The source position in the ECS is s v = (α, 0, d ), where α = s • e e x . See Figure 45. The orientation of the virtual detector in the ECS is characterized by the Euler angles θ = φ = η = 0, i. e. the rotation matrix is the Identity matrix I 3 . Finally, the (virtual) focal distance is d , the vertical offset of the detector is zero and the horizontal offset is α. We can build the projection matrix P v from this geometrical description.

P v = K v R v I 3 | -s v ( 230
)
with

K v =     -d 0 α 0 -d 0 0 0 1     , (231)
R v = I 3 , (232)
s v =     α 0 d     . (233)
We also let

A v = K v R v so that P v = A v | -A v s v .
Note that s p = s v . We keep the subscripts for consistency of notation but both are the triple of coordinates of the source in the ECS.

S T E P 2 :

T H E H O M O G R A P H I C C H A N G E O F V A R I A B L E S .
To derive the change of variables, we start from one point (u, v) and denote r the ray joining the point (u, v) and the source. We search for an expression of (u, v) in terms of the coordinates (U ,V ) of the point of intersection of r with the physical detector. See Figure 50. In homogeneous coordinates, the point is (u, v, 1) (see Appendix 1.3). Any 3D point M in ECS coordinates on the ray r (joining the pixel (u, v) to the source s v ) can be written (U ,V )

M = s v + κA -1 v     u v 1     , κ ∈ R. (234)
Figure 50: One projection and the virtual detector. The orientation of the virtual detector is determined by the epipolar geometry, characterized by the baseline b (red). In Chapter II-2, FBCC were computed from projection data re-sampled on the virtual detector (red), with coordinates (u, v).

In order to implement FBCC directly in the physical detector (blue), with (U ,V ) coordinates, we derive the change of variables from (U ,V ) to (u, v).

This relation is checked by re-projecting M onto the virtual detector:

P v M 1 = A v | -A v s v        s v + κA -1 v     u v 1     1        (235) = A v s v + κ     u v 1     -A v s v (236) = κ     u v 1         u v 1     . ( 237
)
The last line confirms that M projects onto (u, v) on the virtual detector.

Let us assume for now that the same ray r intersects the physical detector in coordinates (U ,V ) (with homogeneous coordinates (U ,V , 1)). We must have

M = s p + κ A -1 p     U V 1     , κ ∈ R. ( 238
)
By identification (and recalling that s v = s p ), we have

    u v 1     A v A -1 p     U V 1     = H     U V 1     , ( 239
)
where we defined the homography H = A v A -1 p . Note that the previous equality is only valid up to a multiplicative constant since it is an equality of homogeneous coordinates. The homography H establishes a correspondence between finite points (U ,V , 1) of the physical detector with finite points (u, v, 1) of the virtual detector. It is more precisely given by

H = A v A -1 p (240) =     -d 0 α 0 -d 0 0 0 1         -e e x - -e e y - -e e z -         | | | e o x e o y e o z | | |          -1 f 0 u 0 f 0 -1 f v 0 f 0 0 1      . ( 241
)
We define some further notation related to H . The (i , j )-coefficient of H is H i j . For j = 1, 2, 3, H • j is the (column) vector H 1 j H 2 j H 3 j T . For i = 1, 2, 3, H i • is the (column) vector H i 1 H i 2 H i 3 T . In the following, we restrict our attention to the finite points (U ,V ) on the physical detector which we note

X =     U V 1     . (242)
The scalar product of X with the i -th line of H is denoted

H T i • X = H i 1 U + H i 2V + H i 3 . We have H X =     H T 1• X H T 2• X H T 3• X         u v 1     . ( 243
)
To express G f in terms of physical coordinates (U ,V ), we would like to express u and v in terms of U and V . Provided that H T 3• X = 0, the above equation provides the desired change of variables

u = H T 1• X H T 3• X , ( 244
) v = H T 2• X H T 3• X . ( 245
)
Before proceeding to the actual change of variables, we need to check if the condition H T 3• X = 0 is compatible with our assumptions. To this end, let Σ denotes the set of points X such that H T 3• X = 0. This condition means that the point (U ,V ) is mapped to a point at infinity in the virtual detector, i. e. the ray from the source to (U ,V ) is parallel to the virtual detector. The set Σ is the intersection of the physical detector and the plane parallel to the virtual detector and containing the baseline b (we call this plane the baseplane). The epipole ε belongs to Σ. But we assumed that the baseline b does not intersect the support of the object so that the epipole is outside the support of the projection. In the next paragraph, we implement the change of variables in the integral expression G f (θ) of Equation 214(subscript i is dropped), which we recall here:

G f (θ) = δ(v -v θ ) g λ (u, v) (u -α) 2 + v 2 + d 2 du dv. (246)
S T E P 3 :

T H E C H A N G E O F V A R I A B L E S .
We prove in Appendix 2.2 that the jacobian of the change of variables of Equations 244 and 245 is

du dv = |Jac| dU dV = | det H | |H T 3• X | 3 dU dV . ( 247
)
By definition of H , we have det H = d 2 / f 2 . For further simplification, we now prove the Lemma 5.

H T 3• X = cos β v cos β p , ( 248
)
where β v and β p are the incidence angles of the ray from s to X with respect to the virtual and physical detector respectively 34 . 34 The notation β p and β v for the incidence angles is only a short-cut. It depends on the detector's variables (U ,V ) and (u, v) respectively.

Proof. From Equation 241, we have

H T 3• X = e e z •          | | | e o x e o y e o z | | |          -1 f 0 u 0 f 0 -1 f v 0 f 0 0 1          U V 1          (249) = e e z •      R T      -U -u 0 f -V -v 0 f 1           (250) = - 1 f e e z •     R T     U -u 0 V -v 0 -f         . (251) (252)
Let us denote α the unit direction of the ray from the source position to the pixel (U ,V ). Since the source has coordinates (u 0 , v 0 , f ) in the 3DCS (O d , u, v , w ), the vector from the source to the point (U ,V ) is l α where

l = (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 . ( 253
)
Its coordinates in the 3DCS (and in the OCS) are (Uu 0 ,Vv 0 ,f ). Left multiplying these coordinates by R T , which is the rotation matrix from OCS coordinates to PCS coordinates, yields the PCS coordinates so that

- 1 f e e z •     R T     U -u 0 V -v 0 -f         = - 1 f e e z • l α (254) = l f cos -e e z , α . ( 255
)
The proof is complete since

• the angle made by the vectorse e z and α is the incidence angle with respect to the virtual detector, namely β v ,

• f /l is the cosine of the incidence angle with respect to the physical detector, namely cos β p .

For geometrical reasons, we note that both cosines are positive.

For reasons that will become clear later, we split the jacobian in two separate pieces

du dv = |Jac| dU dV = | det H | |H T 3• X | 2 | cos β p | | cos β v | dU dV . ( 256
)
Before plugging the new variables in the integral, we note that the weighting factor in Equation 246satisfies

1 (u -α) 2 + v 2 + d 2 = cos β v d . ( 257
)
Symmetrically on the physical detector

cos β p = f (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 . ( 258
)
After some rearrangement, we obtain du dv

(u -α i ) 2 + v 2 + d 2 = | det H | |H T 3• X | 2 f (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 dU dV d ( 259
)
and finally substitute (u, v) variables with their expression in terms of (U ,V ) coordinates in the expression of G f :

G f (θ) = δ H T 2• X H T 3• X -v θ f g λ (U ,V ) (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 | det H | |H T 3• X | 2 dU dV d . ( 260
)
Note that g λ (u, v) is replaced by g λ (U ,V ) because g λ (u, v) refers to the projection image on the virtual detector whereas g λ (U ,V ) refers to the projection image on the physical detector. As announced, the cosine of the incidence angle with respect to the virtual detector is now substituted with the cosine of the incidence angle with respect to the physical detector.

m(U ,V ) = H T 2• X H T 3• X -v θ . ( 261
)
The function m(U ,V ) is smooth for all (U ,V ) on the physical detector (recall that we assumed that the set Σ of singular points does not intersect the physical detector). We use the following theorem (which is presented in the Appendix 1.1, Theorem 23 in a slightly more general version) 35 . 35 This theorem is not a contribution of this work.

Theorem 19 (Composition of δ with a smooth function.). Let m be a smooth function on R 2 and φ a test function. Let Γ = x : m( x) = 0 . We assume that ∇m( x) = 0 for all x ∈ Γ.Then:

R 2 δ (m( x)) φ( x) d x = Γ φ( x) ∇m( x) dσ, ( 262
)
where the integral must be understood in the sense of curve integral.

In view of applying Theorem 19 to our expression, we prove in Appendix 2.3 that

∇m (U ,V ) = (∆ 13 U -∆ 11 ) 2 + (∆ 13 V -∆ 12 ) 2 H T 3• X 2 , ( 263
)
where ∆ i j is the cofactor of the entry (i , j ) of the matrix H . E. g.

∆
H -1 =     -f 0 u 0 0 -f v 0 0 0 1         -e o x - -e o y - -e o z -         | | | e e x e e y e e z | | |         -1 d 0 -α d 0 -1 d 0 0 0 1     , (264)
so that

∆ 13 = (det H ) H -1 31 (265) = (det H ) e o z • - e e x d (266) = - d f 2 e o z • e e x ( 267
)
which vanishes if and only if the baseline is parallel to the physical detector. We excluded this eventuality so ∆ 13 = 0. Equation 263becomes

∇m (U ,V ) = U - ∆ 11 ∆ 13 2 + V - ∆ 12 ∆ 13 2 |∆ 13 | H T 3• X 2 . ( 268
)
The geometrical meaning of the offsets in the square root expression is elucidated in the following Lemma 6 (Epipole coordinates). The coordinates of the epipole ε in the DCS

are (U ε ,V ε ) = ∆ 11 ∆ 13 , ∆ 12 ∆ 13 .
Proof. The ray from the epipole through the source (the baseline) is parallel to the virtual detector in the direction of e e x . In terms of homogeneous coordinates, this means that the finite point U ε ,V ε , 1

T in the physical detector is mapped to the point at infinity 1, 0, 0 T in the virtual detector. In matrix form, this is written

H     U ε V ε 1         1 0 0     , or     U ε V ε 1     H -1     1 0 0     . ( 269
)
The inverse of H is given by

H -1 = 1 det H     ∆ 11 ∆ 12 ∆ 13 ∆ 21 ∆ 22 ∆ 23 ∆ 31 ∆ 32 ∆ 33     T (270) so that (recall that ∆ 13 = 0)     U ε V ε 1     1 det H     ∆ 11 ∆ 12 ∆ 13         ∆ 11 /∆ 13 ∆ 12 /∆ 13 1     , (271)
which ends the proof of the lemma.

Together with Lemma 6, we can apply Theorem 19 to the integral in Equation 260. The |H T 3• X | 2 terms cancel out and we finally obtain

G f (θ) = {m(U ,V )=0} f g λ (U ,V ) (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 | det H | d |∆ 13 | dσ (U -U ε ) 2 + (V -V ε ) 2 (272)
We are two more rearrangements away from the announced result:

• We shall prove that the integration set {m(U ,V ) = 0} is precisely the epipolar line l θ . Since we assumed that H T 3• X = 0 for all (U ,V ) in the physical detector, we have

m(U ,V ) = 0 ⇔ H T 2• X -v θ H T 3• X = 0 (273) ⇔ 0 1 -v θ T H     U V 1     = 0 (274) ⇔ l T X = 0, (275)
where we put l = H T (0, 1, -v θ ) T . As explained in Appendix 1.3, Proposition 27, the line H T (0, 1, -v θ ) T is the image of the line v = v θ through the homography H -1 , i. e. the epipolar line l θ .

• From Equation 267 The final result is

G f (θ) = 1 cos e o z , e e x l θ f g λ (U ,V ) (U -u 0 ) 2 + (V -v 0 ) 2 + f 2 dσ (U -U E ) 2 + (V -V E ) 2 (276)
which ends the proof of Theorem 18.

The integral in Equation 276 integrates the (physical) projection g λ weighted by the cosine of the (physical) incidence angle and by the inverse of the distance to the epipole, along the epipolar line l θ . By parametrizing the 2D unit vectors on the physical detector plane with their angle γ ∈ 0, 2π with the U -axis, we have α(γ) = (cos γ, sin γ) in the DCS. We define γ θ so that α(γ θ ) is the unit vector in the direction of l θ . Equation 276 can be rewritten in the more compact form

G f (θ) = 1 cos e o z , e e x ∞ 0 gλ ε + t α(γ θ ) t dt , ( 277
)
where gλ denotes the cosine-weighted projection. This form of G f will be useful for the comparison with GCC proposed in the next section.

G g (λ, β) = S 2 δ α • β gλ α d α, ( 278
)
whose calculation consists in the following four steps:

• Weight the projection with the cosine of the incidence angle on the detector,

• Compute the 2D Radon transform of the weighted image in the u-direction.

• Differentiate in the v-direction.

• Post-weight the result with the squared cosine.

The original paper of Grangeat focussed on reconstruction. Projections were considered independently from each other. When considering a infinitesimal change in the plane E θ , Grangeat considered planes whose intersections with the detector were parallel (this corresponds to an epipole at infinity) and opted for a differentiation in the orthogonal v direction. When considering a pair of projections, the epipoles become finite points in the detectors' planes. An infinitesimal change in the plane E θ results in a epipolar line on the detectors. It seems more natural to differentiate with respect to the angle of the epipolar line (see Figure 51 right) rather than the v-direction. The next proposition addresses this idea.

The same set-up applies: We consider an epipolar geometry, characterized by the baseline is b and an ECS. We restrict our attention to one projection g λ . An epipolar plane E θ intersects the detector in an epipolar line l θ . The normal to E θ is denoted β(θ). The cosine-weighted projection is denoted gλ . As above, unit vectors in the detector plane are denoted α(γ) where γ is the polar angle with respect to the u-axis of the DCS. With this notation, α(γ θ ) denotes the unit direction of the epipolar line l θ . We introduce the shorter notation G g (θ) = G g λ, β(θ) . We can now state the following Proposition 20 (Epipolar Grangeat theorem). For all θ

G g (θ) = 1 cos 2 θ ∂ ∂γ ∞ 0 g ε + t α(γ) t dt γ=γ θ (279)
Compared to the previous 4-step procedure, we note the following changes:

• An extra weighting by the distance to the epipole must be added. • Instead of computing the Radon transform, a fanbeam parametrization (centred in ε) of the line l θ is used,

• The differentiation step is now with respect to the epipolar angle (i. e. the angle of the fanbeam parametrization).

Proof. We start from Equation 77in Section 1.7 (with (U ,V ) coordinates), replace U u +V v with x for ease of notation, and note that

s λ • β(θ) = ε• β(θ).
We obtain

G g (θ) = δ (U u + V v -s) • β(θ) g (U ,V ) dU dV (280) = δ x • β(θ) -s λ • β(θ) g ( x) d x, ( 281
) = δ x • β(θ) -ε • β(θ) g ( x) d x. ( 282
)
We substitute x with polar coordinates centred in the epipole ε, i. e.

x = ε + t α(γ) t ≥ 0, γ ∈ 0, 2π . (283)
The Jacobian determinant of this change of variables is t . Using the homogeneity property of δ , we obtain

G g (θ) =
β D • α(γ) = ± β D sin(γ -γ θ ). ( 287
)
We use the homogeneity of δ and apply Proposition 24 of Appendix 1.1 to obtain

G g (θ) = 1 β D 2 ∂ ∂γ +∞ 0 gλ ε + t α(γ) t dt
) 291
The proof is complete.

The Equation 279 is simply another way of evaluating G g (β). The output is exactly the same as with Equation 67. In fact, there are as many ways of computing G g (β) as points on the epipolar line l θ . Each point on l θ can be considered as the epipole of some epipolar pair, and the same approach applies.

One of the objectives of this chapter was to translate the intrinsic relation between FBCC and GCC of Proposition 16 in the physical detector. In Proposition 16, G f was differentiated with respect to the angle θ of the epipolar plane. In the expression of G g in the physical detector of Proposition 20, the derivative is taken with respect to the angle γ θ made by the epipolar line in the physical detector. These two angles are obviously connected to each other. In the next section, we investigate this connection and complete this study.

T H E I N T R I N S I C R E L AT I O N B E T W E E N F B C C A N D G C C O N T H E P H Y S I C A L D E T E C T O R

We recall the Proposition 16 which relates G f and G g in spherical coordinates (with slight change in notation). For all θ ∈ -π 2 , π 2 ,

G g θ = ∂ ∂θ G f (θ) θ= θ. ( 292
)
This relation must be true if the expression of G g and G f are replaced with their equivalent expressions on the physical detector (Proposition 20 and Equation 277 respectively). We start from the RHS of Equation 292, use the expression of G f (θ) in the physical detector (Equation 277), and apply the chain rule:

∂ ∂θ G f (θ) θ= θ = ∂γ θ ∂θ ∂ ∂γ 1 | cos( e o z , e e x )| ∞ 0 g ε + t α(γ) t dt γ=γ θ . ( 293
)
We then have the following Lemma 7.

∂γ θ ∂θ = cos( e o z , e e x ) cos 2 θ . ( 294
)
The proof can be found in Appendix 2.4. Finally, Equation 293becomes

∂ ∂θ G f (θ) = ± 1 cos 2 θ ∂ ∂γ ∞ 0 g ε + t α(γ) t dt γ=γ θ , (295)
which is precisely (up to the sign) the expression of G g in the physical detector, derived in Proposition 20. The sign is related to the individual projection from the pair. The cosine term cos( e o z , e e x ) will be positive for one and negative for the other, hence the sign flip.

N U M E R I C A L S I M U L AT I O N

This section presents preliminary results, which validate numerically the theoretical results presented above. We simulated one pair of projections. The angular gap between the projections was set to 30 • . The detector was placed at 200 mm in both cases. It was perfectly aligned in the first projection and tilted around its v-axis by 3 • in the second projection. The size of the detector was 256mm × 256mm with a pixel size of 0.25 mm (512 × 512 pixels). The resulting half cone-angle was 18 • approximately.

We simulated the projections of a thorax Forbild phantom, which was placed at the origin O of the PCS, at equal distance from the source and the detector. The two projection images are presented in Figure 54. They were computed analytically from the definition of the phantom.

For each projections i = 1, 2, we computed

• G f ,i (θ) with Equation 277, • G g ,i (θ) with Equation 279,

• The derivative of G f with respect to the angle of the epipolar plane θ.

The two analytic noise free projections, along with their respective projection matrices are consistent, so that the functions G f ,1 (θ) and G f ,2 (θ) must be equal. Discrepancy between the two signals would necessarily stem from the implementation.

N U M E R I C A L V A L I D AT I O N O F T H E F B C C I N T H E P H Y S I C A L D E T E C - T O R
We sampled the function G f ,i (θ) (i = 1, 2), for a regular sampling of θ ∈ [-36 • , 36 • ] (the extreme values were manually chosen so that the whole support of the projection images was scanned). The Figure 55 presents the plots. Both curves match perfectly, which confirms the expression of G f in the physical detectors. For comparison purposes, we also present in Figure 55 the same plots with G g ,i .

In order to evaluate the robustness of the implementation to noise, we degraded the projections with Poisson noise, according to the following procedure. The number of photons received by a pixel when no object is placed was set to I 0 = 10 4 . The densities of the phantom were multiplied by the attenuation coefficient of water (0.01879mm -1 at 75 keV). Then a Poisson distribution was drawn from the computed intensities. The results are presented in Figure 56. Both curves still strongly agree in the case of FBCC. The matching is less tight in the case of GCC. The differentiation step in G g is more sensitive to noise. This explains the higher discrepancies.

T H E R E L AT I O N B E T W E E N F B C C A N D G C C
In Figure 57, we illustrate the relation

G g (θ) = ∂ ∂θ G f (θ) (296)
by plotting G g (θ) on one hand and ∂/∂θG f (θ) (according to Equation 293) on the other hand, computed on one projection of the pair (in this case, the first one). C

O N C L U S I O N A N D F U T U R E W O R K S U M M A R Y O F A C H I E V E M E N T S
This work has achieved some significant results, both theoretical and practical. On the theoretical side, we introduced in Chapter II-1 new consistency conditions for circular conebeam projections. These conditions were derived by revisiting the Radon intermediate functions from the perspective of consistency. We also focussed on pre-existing fanbeam consistency conditions and used them in a pairwise cone-beam context. This has lead us to a thorough comparison of these DCC and previously published Grangeat based consistency conditions. On the practical side, we applied the fanbeam DCC to the geometric calibration problem and proposed a fully on-line calibration of a micro-CT system.

F U T U R E W O R K
The end of a PhD period is quite frustrating because one has to spend several months writing a manuscript and in doing so, realizes how much is left to be done with no time for doing it! Here are a few research directions which may be worth investigating. First, on the DCC side, general necessary and sufficient conditions for cone-beam data are still defying the CT community. I keep thinking that a geometrical approach as the one presented in [START_REF] Clarkson | Symmetry properties of an imaging system and consistency conditions in image space[END_REF] and summarized in Section I-2 could be promising. Also, necessary and sufficient conditions exist for cone-beam projections with the extra assumption that the support of the object does not intersect the plane of the trajectory of the source. This kind of assumption is typically to avoid division by zero, which the projective geometry can handle. Revisiting the conditions with the projective geometry toolbox may prove to be effective.

In terms of applications, and specifically of geometric calibration, this work proposed a proof-of-concept but the fully on-line calibration of a Carm is still on-going research. The extension of the proposed method to projection-specific calibration is, in theory, possible. The resulting high dimension of the optimization problem may, in practice, be difficult to handle. A first step would consist in geometrically registering a pair of projections (i. e. determining the geometry of one projection with respect to the other, based on pair-wise consistency). Finally, the field of computer vision has proposed numerous geometric calibration methods, including self-calibration of the intrinsic matrix. These methods are usually based on the analysis of sharp points in the images. Natural images essentially differ from X-ray projections images because the line integral forward model is regularizing the object function. Nevertheless, a human eye can see sharp edges in a radiograph. The design of hybrid calibration methods, mixing computer vision algorithms with consistency analysis may be worth investigating.

Part III

A P P E N D I X A A P P E N D I X A : A M AT H E M AT I C A L T O O L B O X 1.1 A S H O R T I N T R O D U C T I O N T O D I S T R I B U T I O N S
In this section, we give a brief overview of distributions. Exhaustiveness is out of reach in the appendix of a PhD dissertation. Nevertheless, it is necessary for the sake of completion, to introduce the properties of very common distributions, which are extensively used in this manuscript. These distributions are the Dirac delta distribution δ, its derivative δ , the Cauchy principal value of 1/s and the Hadamard's finite part of 1/s n . These distributions play a central role in tomography and many of the integrals involved are only defined in the sense of distributions (though this is not always clearly mentioned in the literature or in this manuscript). The theory of distribution has this particularity that everything goes as usual while the object being studied are simply not defined in the usual world. When we defined the ramp filter (as the inverse Fourier transform of the function |σ|) or the Hilbert filter (as the inverse Fourier transform of -i sign(σ)) and carry out complex calculations as usual, we omitted that in the usual world (i. e. the world of functions), the functions |σ| and sign(σ) are not integrable and have no inverse Fourier transform. The theory of distributions was precisely designed by Laurent Schwartz 36 to mathematically guaran- tee that all these formal derivations remain valid. The DCC derived in the first contribution are essentially a generalization of the ramp filter (seen as the derivative of the Hilbert kernel) to higher order of differentiation. Finally, the third contribution relies crucially on a theorem which rigorously defines the composition of the δ distribution with smooth functions. This introduction to distributions was written with the help of various sources of information, mostly [START_REF] Strichartz | A Guide to Distribution Theory and Fourier Transforms[END_REF], lecture notes from the web, wikipedia, etc...

The space of test functions D

We focus our presentation on the one-dimensional case. We start with the following Definition 7. A test function is a C ∞ function defined on R whose support 37 37 The support of a function φ is the smallest closed set containing all x ∈ R such that φ(x) = 0. Applying any translation or dilation to φ yields infinitely many other examples of test functions. It is very easy to prove that D is a vector space. One finally wants to define a topology on this vector space. This point is the key of the whole theory and is a bit difficult. We limit our presentation to the characterization of the sequential convergence in D.

Definition 8 (Sequential convergence in D). For φ k k∈N , φ ∈ D, φ k D -→ φ if there exists a compact set K ⊂ R such that Supp φ k ⊂ K ∀k ∈ N, (298)
∂ m ∂x m φ k -φ ∞ k -→ 0 for all m ∈ N. (299)
The algebraic dual space of D is the set of all linear forms on D. Distributions are the linear forms on D which are continuous for the topology induced by the convergence defined above. It is well known that the continuity of a linear form is equivalent to the local continuity in 0 and that the continuity is equivalent to the sequential continuity. We then have the following characterization of continuous linear forms. We say that the concept of distribution is an extension of the concept of usual functions because usual functions can be identified with elements of D . More precisely, we denote L loc 1 the set of locally integrable functions, i. e. integrable on every compact set 38 . For every f ∈ L loc 1 , we define T f on 38 In this space, we identify two functions which agree almost everywhere.

Proposition 21. Let T be a linear form on D. T is continuous if and only if

φ k D -→ 0 ⇒ T (φ k ) → 0. ( 300
D by

∀φ ∈ D, T f φ = f (x)φ(x) dx. ( 301
)
The mapping T f is clearly linear. To prove the continuity, we consider a sequence φ k D -→ 0 and the compact K as in Definition 8. We have

|T f (φ k )| = K f (x)φ k (x) dx ≤ φ k ∞ K | f (x)| dx. ( 302
)
In the last expression, φ k ∞ → 0 by assumption and the integral is well defined because f is locally integrable. Finally |T f (φ k )| → 0. Consider now two functions f and g ∈ L loc 1 such that T f = T g i. e. ∀φ ∈ D, T f (φ) = T g (φ). We then have

( f -g )(x)φ(x) dx = 0 ∀φ ∈ D, (303)
so that f = g almost everywhere, i. e. f = g in L loc 1 . The mapping from L loc 1 to D is injective. The distributions T f stemming from L loc 1 functions are called regular distributions. We give important examples of such regular distributions.

• T H where H is the Heavyside function i. e. the indicator function of R + .

• T sign where sign is the signum function i. e. sign(x) = H (x) -H (-x).

• T 1 where 1 is the constant function, which is one (almost) everywhere.

Of course, there are distributions which are not regular (they are called singular). We now give some examples of such singular distributions and focus on those of interest in this work. In the following T (φ) will be denoted < T , φ > to fit the conventional distribution notation.

T H E D I R A C D I S T R I B U T I O N

The first example is the Dirac δ distribution, defined by < δ, φ >= φ(0). It cannot be a regular distribution since as a function, it would be zero everywhere but in zero, so it would be the zero function in L loc 1 and would satisfy δ(x) dx = 1 (take any φ ∈ D such that φ(0) = 1), which is impossible. We will say more about the Dirac distribution in Section 1.1.5.

P R I N C I P

A L V A L U E O F 1/x : pv(1/x )
The second example is what we called the Hilbert kernel. The function 1/x has a singularity in zero so is not locally integrable. It can not be a regular distribution. Nevertheless, one would desire to have an equivalent to the function 1/x. That is the purpose of the distribution pv 1

x . Proposition 22 (The principal value). For all φ ∈ D, the map

< pv 1 x , φ >= lim →0 + - -∞ φ(x) x dx + +∞ φ(x) x dx . ( 304
)
is a distribution. Proof. Provided, the limit exists, it is clear the the above equation defines a linear map. Since the support of φ is compact, the integral in ±∞ is not a problem. We write

- -∞ φ(x) x dx + +∞ φ(x) x dx = +∞ φ(x) -φ(-x) x dx (305)
The integrand in the RHS integral can be continuously extended in 0 since its limit in zero is 2φ (0). Hence the existence of the limit. To prove the continuity, we consider a sequence (φ k ) and a compact K like in Definition 8.

From the Taylor expansion of φ, we have

φ(x)-φ(-x) x = 2φ ( x) for some x. Thus +∞ 0 φ k (x) -φ k (-x) x dx ≤ φ k ∞ K 1 dx k ----→ 0 (306) F I N I T E P A R T O F 1/x n ,
(x) = φ(0) + xφ (0) + x 2 φ (0)/2 + .... ( 307
) Let A > 0 such that Suppφ ⊂ [-A, A] and let K = [-A, -] [ , A]. We have K φ(x) x 2 dx = K φ(0) x 2 dx + K φ (0) x dx + K φ (0) 2 dx + ..., (308)
The terms of order 2 and higher are no problem since they integrate a polynomial over a bounded set. The term of order 1 is zero (it is the integral of an odd function on a symmetric domain). The term of order zero is K φ(0)

x 2 dx = - 2φ(0) A + 2 φ(0) . ( 309
)
The term 2 φ(0) is obviously diverging. So we define the finite part distribution by subtracting the divergent part. More precisely, we have the Definition 9 (Finite part of 1/x 2 : fp(1/x 2 )).

< fp 1

x 2 , φ >= lim →0 |x|> φ(x) x 2 dx -2 φ(0) . ( 310
)
We define similarly the distribution fp(1/x n ). For instance, Definition 10 (Finite part of 1/x 3 : fp(1/x 3 )).

< fp 1

x 3 , φ >= lim →0 |x|> φ(x) x 3 dx -2 φ (0) . (311)

Derivative of a distribution

One of the nice features of distributions is that we can define a differentiation which is not restricted to differentiable functions. All distributions have a derivative, which itself is a distribution, so that all distributions are infinitely differentiable. The definition given above coincides with the usual differentiation for the differentiable function. Let f be a differentiable function. A simple integration by parts gives

k ⊂ Suppφ k ⊂ K . Finally < DT , φ k > = < T , φ k > k ----→ 0. ( 313
) T DT DT H = δ DT sign = 2δ DT 1 = 0 < Dδ, φ > = -φ (0) < D n δ, φ > = (-1) n φ (n) (0) < Dpv 1 x , φ > = -fp 1 x 2 < Dfp 1 x 2 , φ > = -fp 1 x 3
< T f , φ > = f (x)φ(x) dx (314) = f (x)φ(x) +∞ -∞ -f (x)φ (x) dx. ( 315
)
The first term vanishes due to the compact support of φ. Hence, the distributional derivative of a regular distribution is the regular distribution associated to the derivative of the function:

DT f = T f .
We give the derivatives of usual distributions in Table 5.

Tempered distributions and Fourier transform

We would like to define the Fourier transform of a distribution as we defined the differentiation, i. e. for φ ∈ D

< F T , φ >=< T , F φ > . (316)
Unfortunately, the Fourier transform of a compactly supported C ∞ function is not compactly supported, so that the RHS of the above equation does not make sense. We need to extend the space of test functions to a space which is stable under the Fourier transform. This is precisely the interesting feature of the Schwartz space S (R). The Fourier transform is a bijective bi-continuous mapping of S (R). We call tempered distribution, a continuous linear form on S . The vector space of tempered distributions is denoted S . Since D ⊂ S , we have S ⊂ D . By extending the space of test functions, we can now define the Fourier transform. Definition 12. Let T ∈ S be a tempered distribution. The Fourier transform of T is given, for all φ ∈ S ,

< F T , φ >=< T , F φ > . ( 317
)
Once again, the definition is compatible with the standard Fourier transform on functions. Let T f be a regular distribution, with f ∈ L 1 .

< F T f , φ > =< T f , F φ > (318) = f (x) φ(y) e -2i πx y dy dx (319) = φ(y) f (x) e -2i πx y dx dy (320) =< T F f , φ >, (321)
where the inversion of signs is an application of Fubini's theorem, legal because

f (x)φ(y) dx dy = f (x) dx φ(y) dy < ∞ (322)
Before giving the Fourier transform of some usual distributions, we formally derive, as an illustrative example, the Fourier transform of the Dirac distribution.

< F δ, φ > =< δ, F φ >, (323)
= F φ (0), (324)
= φ(x) dx, ( 325
) =< 1, φ > . (326)
Finally, F δ = 1. The Fourier transforms of some common distributions are given in Table 6. No proof is given, due to some technicalities which would make the presentation too long.

We complete this section with the Fourier transform of δ and fp 1 x 2 . To this end, it is necessary to say of few words about multiplication of distributions. Strictly speaking, the multiplication of two distributions cannot be defined. We can though define the multiplication of a distribution T by a smooth function ϕ by, again, transferring the operation from the distribution to the test function. More precisely where ϕφ ∈ D. In particular, one can define the multiplication by the identity function x:

< ϕT , φ >=< T , ϕφ >, (327)
T F T F δ = 1 F T H = δ 2 -i 2π pv 1 x F T sign = i π pv 1 x F T 1 = δ F pv 1 x = -i πsign(x)
< xT , φ >=< T , xφ > . ( 328
)
Remembering that the derivative of the Fourier transform (F φ) (x) is -2i π times the Fourier transform of xφ(x), we derive the Fourier transform of δ as

< F δ , φ > =< δ , F φ > (329) = -< δ, F φ > (330) = 2i π < δ, F (xφ(x)) > (331) =< F δ, 2i πxφ > (332) =< 2i πx, φ > . ( 333
)
Similarly, we derive the Fourier transform of fp 1 x 2 :

< F fp 1 x 2 , φ > =< fp 1 x 2 , F φ > (334) = -< Dpv 1 x , F φ > (335) =< pv 1 x , F φ > (336) = -2i π < pv 1 x , F (xφ(x)) > (337) = -2i π < F pv 1 x , xφ > (338) = -2i π < -i πsign(x), xφ > (339) = -< 2π 2 |x|, φ > . ( 340
)
Finally, we can now gather the main results which justify the definitions of the four distributions introduced in Chapter I-1: δ, δ , the Hilbert kernel h and the ramp filter ρ. They were simply defined as the inverse Fourier transform of some functions which ... had no inverse Fourier transform, in the usual sense. With the distribution theory, we may rigorously write δ(x) = e 2i πxσ dσ (341) The Dirac distribution was defined in the previous section, and its most important properties (derivative, Fourier transform, etc.) were derived. Nevertheless, we give two more important results which are needed in Chapter II-3. They are both related to the composition of δ with smooth maps.

δ (x) = 2i πσ e 2i πxσ dσ (342) h(x) = (-i sign(x)) e 2i πxσ dσ (
By a simple change of variable in the definition above of δ, one can easily prove that, for all a > 0

δ(ax) = δ(x) a . ( 345
)
This must be understood in the sense of distribution. For all φ ∈ D δ(ax)φ(x) dx = φ(0)

a . ( 346
)
The composition of δ can be generalized to more complex C ∞ (R n ) functions, with the following theorem, whose proof is not presented (see [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I[END_REF]) Theorem 23 (Composition of the Dirac δ distribution with smooth function.). Let m be a smooth function on R n and φ a test function. Let Γ = x : m( x) = 0 . We assume that ∇m( x) = 0 for all x ∈ Γ.Then:

R n δ (m( x)) φ( x) d x = Γ φ( x) ∇m( x) dσ, ( 347
)
where dσ is the Euclidean surface measure on Γ.

Of course, Equation 345 is a particular case of the theorem, with m(x) = ax. The proof of this theorem is difficult and out of scope of this Appendix. We prove the theorem in the particular case of a one-dimensional function m which has only one root in zero (i. e. m(x) = 0 ⇔ x = 0). We also assume this root is simple: m (0) = 0. By assumption, we first note that the function τ(x) = m(x)/x, x = 0 can be continuously extended in 0 with τ(0) = m (0). Writing δ as the inverse Fourier transform of 1 and changing the Fourier variable, one obtains

which is conform to intuition.

I N T E G R AT I O N O V E R S P H E R E S

Though there is a general theory on the integration of a function over unit spheres, we only give the formula for S 1 and S 2 . The goal is to compute the integral of a function f defined on R (resp. R 2 ) on the sphere S 1 (resp. S 2 ).

Proposition 25 (Integration over S 1 ). Proof. This is an immediate consequence of the general formula for surface integral:

S f ( x) d x = f ( x(x 1 , x 2 )) ∂ x ∂x 1 × ∂ x ∂x 2 dx 1 dx 2 , (363)
where x(x 1 , x 2 ) is a parametrization of the surface S 2 , × denotes the crossproduct of two 3D-vectors and • is the standard 2-norm on R 3 . Here, the parametrizations of the sphere is given by: α(θ, ϕ) = (cos θ sin ϕ, sin θ sin ϕ, cos ϕ), (364) α(θ, ϕ) = (cos θ cos ϕ, sin θ cos ϕ, sin ϕ), (365)

for the latitude-longitude and colatitude-longitude cases respectively. The computation of the surface element ∂ x ∂x 1 × ∂ x ∂x 2 is straight-forward.

T H E B A S I C S O F P R O J E C T I V E G E O M E T R Y

In the CT community, projection matrices are very commonly used. But the underlying concepts of the projective geometry seem to be less familiar (at least, few research articles use these concepts). When considering nice acquisition geometries (circular orbit of the source, detector perfectly aligned), elementary geometry may perfectly do the job and projective geometry is not necessary. If the trajectory and the respective position of the detector and the source are more general, computations with elementary geometry become intractable. The concise shape of linear algebra tools to describe one projection (the projection matrix) or transformation between a pair of projections is worth the journey in the projective geometry.

Projective geometry is the geometry that best describes human vision. Heuristically, the human eye which looks at railways sees lines converging at the horizon. Projective geometry gives a rigorous sense to such statements as "parallel lines meet at infinity". This introduction is concise and focuses on the concepts needed in the course of the text. We give a formal definition of projective space and a geometric interpretation, in relation with the concept of homogeneous coordinates. Then we introduce homographies. They are the equivalent for projective spaces of non singular linear maps for vector spaces. In particular, we illustrate this with the computation of images of lines through homography in the projective 2D plane.

Projective spaces

The following definition is the true definition of a general (real) projective space.

Definition 13 (Projective space). The real projective space of dimension n ∈ N is the quotient where denotes the equivalence relation (x 0 , ..., x n ) (x 0 , ..., z n ) iif ∃λ = 0 such that (x 0 , ..., x n ) = (λx 0 , ..., λx n ).

P n (R) = R n+1 -0 (366)
A rephrasing of the definition is that a n-dimensional real projective space is the set of all lines in R n+1 through the origin O. Indeed, two points on such a line are equivalent. The problem with this formulation (an element of P n (R) is a line) is that we want to define lines in a projective space just as we do in a affine space. We would then define lines made of... lines. To avoid confusion, we need another interpretation of the definition, closer to the usual geometry. For the sake of clarity, we concentrate on the real projective plane (i. e. P 2 (R)). We start from the 3D space R 3 , endowed with the canonical coordinate system. The situation is illustrated in Figure 58. We let Π 1 and Π 0 denote the affine planes with equation z = 1 and z = 0 respectively. The origin O may be seen as the X-ray source and the plane Π 1 as the plane of the detector. A line through the origin (an X-ray) may or may not intersect Π 1 . If it does, a 3D point on such a line have coordinates (x, y, z) with z = 0 and the line intersects Π 1 in (x/z, y/z, 1). We then have a one-to-one mapping between these lines and the points of the plane Π 1 . We then conclude that the projective plane contains the affine plane Π 1 . We now turn to those lines passing through the origin and not intersecting Π 1 . They are precisely the lines contained in Π 0 . Points in Π 0 have coordinates (x, y, 0). They can be associated to the points of the semi-unit-circle (x/ x 2 + y 2 , y/ x 2 + y 2 , 0) with x > 0 for example.

With the above representation, a triple (x, y, z) ∈ R 3 can be seen as a point of the projective plane. If z = 0, this point is finite and belongs an affine plane. The point is said at infinity if z = 0. There is a one-to-one mapping between points at infinity and half of a unit circle in Π 0 . Finally, we have described the projective plane as an affine plane (finite points) plus a half-circle of points at infinity. This interpretation may be more familiar to a non-mathematician reader. Nevertheless, one should keep in mind that the nature of one particular point (finite or at infinity) is not intrinsic but very much depends on the choice of a coordinate system. E. g. , interchanging the roles of z and x would send the point (0, 1, 1) at infinity.

Points and lines in the projective plane

The equation of a line in an affine space with coordinates (x, y) can always be written

ax + b y + c = 0, ( 367
)
for some l = (a, b, c) ∈ R 3 with a 2 + b 2 = 1 and c ≥ 0. We use l to denote the line. The vector (a, b) is the unit normal to the line and c is its distance to the origin. We now consider the affine points (x, y) as points x of the projective plane x = (x, y, 1), so that Equation 367 can be written in scalar product form (a, b, c) • (x, y, 1) = 0 or l • x = 0.

(368)

Of course, the equation above characterizes the fact that x belongs to the line l if the relation is independent of the representative in the equivalence class of x. In other words, if one has (a, b, c) • (x, y, 1) = 0, then for all λ = 0, one must also have (λa, λb, λc) • (x, y, 1) = 0, which is obviously true. We can now ask if the line l , seen as a line in the projective plane, has a point at infinity. Such point has its third coordinate z = 0. We solve the equation

ax + b y + c × 0 = 0, (369)
up to a multiplicative constant. The set of solutions is of course R(-b, a, 0). This is the projective point at infinity (-b, a, 0), which corresponds to the direction of the line l . The situation is illustrated in Figure 59. This is independent of c, so that all lines parallel to l meet at the same point at infinity. And this is why we say that parallel lines meet at infinity.

I M P O R TA N T R E M A R K
The careful reader may have noticed that lines can be represented in the exact same way as points are, i. e. with a triple of real numbers, up to a multiplicative constant. This remark is the basis of a very fruitful correspondence between points and lines in a projective plane, called duality. The duality generalizes to higher dimensions (lines being hyperplanes in higher dimensions).

Illustration with a projection matrix

Let P = K R I |s be the projection matrix associated to a tomographic device in a 3D Patient coordinate system. s is the position of the X-ray source. By construction, if a 3D point (X , Y , Z ) is projected onto a 2D point (u, v) on the detector, we have the relation

    u v w     = K R I | -s        X Y Z 1        (370)
Let us search for points of the 3D space M = (X , Y , Z ), different from the source position s, which project to points at infinity in the detector projective plane. For such points, we have w = 0. From Chapter I-2, we know that the rows of R are the coordinates of the vectors e o

x , e o y and e o z in the PCS and that K has the form

K =     f 0 u 0 0 f v 0 0 0 1     . (371)
From Equation 370, we have

w = R T 3• (M -s), ( 372
)
which vanishes if and only if 39 Ms is orthogonal to the focal axis (O, e o z ). 39 Remember we excluded the source position s.

Points M which satisfies this condition are precisely the points in the plane parallel to the detector plane and containing the source position. In other words, ray that are parallel to the detector are not detected. All other points in space satisfy w = 0 and project to a point (u, v, w) (u/w, v/w, 1). The situation is depicted in Figure 60.

Plane homography

In this section, we introduce homographies between two projective planes. Mathematically, homographies (or projective applications) are the geometric transforms which preserve the projective structure, i. e. which maps lines to lines. A complete introduction to projective applications is useless here (see [START_REF] Hartley | Multiple View Geometry in Computer Vision Second Edition[END_REF] for a complete treatment from a Computer Vision perspective or [START_REF] Samuel | Projective geometry[END_REF] for an abstract mathematical treatment). We simply define a homography as a transformation which maps homogeneous coordinate triples (see Figure 61). It is represented by a 3 × 3 matrix H so that Note that the matrix is defined up to a multiplicative constant (as one is now familiar with). It is necessary that H is non-singular, for otherwise, some point x ∈ null H would project to H x = 0 which is not a triple of homogeneous coordinates 40 . Re-sampling projection data into a virtual de- tector, as was done in Chapter II-2 provides with a nice archetypal example of homography. This is illustrated in Figure 61. The value of the projection at pixel (U ,V ) should be transferred to the pixel (u, v) of the virtual detector.

    u v w     =   
The following easy proposition is needed in Chapter II-3.

Proposition 27 (Image of a line through an homography). Let H be a 2D homography and l = (a, b, c) T be the representation of a line. Then, the image of l through the homography H is the line l = H -T l .

Proof. By definition, l = H x, ∀x s.t. l T x = 0(i. e. x ∈ l ) . Let H x ∈ l , then:

H -T l T H x = l T H -1 H x = l x = 0 (375)
Conversely, if H -T l T y = 0, then l T H -1 y = 0 and H -1 y ∈ l so that y = H x for some x ∈ l .

so that the Fourier transform f is homogeneous of degree 1 (note that f is not strictly linear but piecewise linear since the constraint is only for a > 0). For all a > 0, one has

By linearity of the inverse Fourier transform, one finally gets that f is a linear combination of δ and ρ. Finally, δ and ρ are clearly linearly independent, which ends the proof.

C O M P U TAT I O N O F T H E J A C O B I A N

In this section, we prove Equation 247

|Jac| = | det H | |H T 3• X | 3 (383)
Proof. From the definition of the change of variables, we have

∂u ∂U = H 11( H T 3• X )-H31(H T 1• X ) (H T 3• X ) 2 ∂v ∂U = H 21( H T 3• X )-H31(H T 2• X ) (H T 3• X ) 2 ∂u ∂V = H 12( H T 3• X )-H32(H T 1• X ) (H T 3• X ) 2 ∂v ∂V = H 22( H T 3• X )-H32(H T 2• X ) (H T 3• X ) 2 . ( 384
)
= 1

H T 3• X 3 det H . ( 390
)

C O M P U TAT I O N O F T H E G R A D I E N T O F φ

In this section, we prove Equation 263∇φ (U ,V ) = (∆ 13 U -∆ 11 ) + (∆ 13 V -∆ 12 )

H T 3• X 2 . ( 391
)
From the definition of φ, we have

∂φ ∂U = H 21 H T 3• X -H 31 H T 2• X H T 3• X 2 (392) ∂φ ∂V = H 22 H T 3• X -H 32 H T 2• X H T 3• X 2 . ( 393
)
Hence

∇φ(U ,V ) 2 = 1 H T 3• X 2 H 2 21 H T 3• X 2 -2H 21 H 31 H T 3• X H T 2• X + H 2 31 H T 2• X 2 + H 2 22 H T 3• X 2 -2H 22 H 32 H T 3• X H T 2• X + H 2 32 H T 2• X 2 . ( 394
)
After expanding all the terms in the denominator and expressing the obtained expression as a polynomial in (U ,V ), we easily obtain:

∇φ(U ,V ) 2 = 1 H T 3• X 4 ∆ 13 U -∆ 11 2 + ∆ 13 V -∆ 12 2 . ( 395
)
Taking the square root of that final equation yields the desired result. In order to derive this relation between the angle θ (angle of the epipolar plane E θ with the plane (O, e e

x , e e z )) and the angle γ (angle of the epipolar line l θ with the U -direction), we write that the line v = v θ in the virtual detector is the image of the line l θ in the physical detector, which in matrix form reads

H T     0 1 -v θ     =     a b c     , (397)
where (a, b, c) is the representation of the line l θ and v θ = d tan θ. In such line representation, (a, b) are the coordinates of a vector orthogonal to the epipolar line l θ . We first prove that Lemma 8. a 2 + b 2 is a constant: for all θ,

a 2 + b 2 = d 2 f 2 . ( 398
)
Proof. Note that d 2 / f 2 = det H for further references. By writing

H T = K -T p R p R 12 K T v , (399)
where β D denotes the orthogonal projection of β(θ) onto the physical detector, i. e. onto the first two components of the OCS e o x and e o y . Taking the squared norm and using Equation 290 and the explanation thereafter, we finally obtain

a 2 + b 2 = d 2 f 2 cos 2 θ cos 2 θ = d 2 f 2 . ( 403
)
The proof is complete.

I- 3 D

 3 ATA C O N S I S T E N C Y C O N D I T I O N S : A R E V I E W 3.1 The range of the Radon transform . . . . . . . . . . . . . . . . 3.2 In divergent beam geometry . . . . . . . . . . . . . . . . . . . I-4 G E O M E T R I C C A L I B R AT I O N O F T O M O G R A P H I C S Y S T E M 4.1 Off-line geometric calibration methods . . . . . . . . . . . . . 4.2 On-line geometric calibration methods . . . . . . . . . . . . .II C O N T R I B U T I O N S II-1C O N E -B E A M D C C F O R C I R C U L A R T R A J E C T O R Y 1.1 New zero-order cone-beam DCC for a circular trajectory . . 1.2 New general cone-beam DCC for a circular trajectory . . . .II-2C O N S I S T E N C Y -B A S E D C A L I B R AT I O N O F A C O N E -B E A M M I C R O -C T 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Notation and background . . . . . . . . . . . . . . . . . . . . . 2.4 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Numerical experiments on simulated data . . . . . . . . . . . 2.6 Numerical experiments on our µ-CT system . . . . . . . . . . 2.7 Conclusion and future work . . . . . . . . . . . . . . . . . . . 2.8 Investigation of the valley-shaped cost function . . . . . . . .II-3F B C C A N D G C C : A C O M P A R AT I V E S T U D Y 3.1 The epipolar geometry of one pair of projections . . . . . . . 3.2 A detector-less version of pair-wise DCC . . . . . . . . . . . . 3.3 Fan-beam DCC in the physical detector . . . . . . . . . . . . 3.4 Grangeat based consistency conditions and the angular form xi xii C O N T E N T S 3.5 The intrinsic relation between FBCC and GCC on the physical detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . C O N C L U S I O N A N D F U T U R E W O R K III A P P E N D I X A A P P E N D I X A : A M AT H E M AT I C A L T O O L B O X 1.1 A short introduction to distributions . . . . . . . . . . . . . . 1.2 Integration over spheres . . . . . . . . . . . . . . . . . . . . . . 1.3 The basics of projective geometry . . . . . . . . . . . . . . . . B A P P E N D I X B : A U X I L I A R Y P R O O F S 2.1 Homogeneous functions of degree -2 . . . . . . . . . . . . . . 2.2 Computation of the Jacobian . . . . . . . . . . . . . . . . . . . 2.3 Computation of the gradient of φ . . . . . . . . . . . . . . . . 2.4 Computation of ∂γ/∂θ. . . . . . . . . . . . . . . . . . . . . . . B I B L I O G R A P H Y A C R O N Y M S PCS Patient Coordinate System. (O, e x , e y , e z ). OCS Oriented Coordinate System. (O, e o x , e o y , e o z ). Given a projection matrix P = K R[I | t ], the OCS is obtained by rotating the PCS by the rotation R. ECS Epipolar Coordinate System. (

1

 1

Figure 1 :

 1 Figure 1: Three X-ray systems.

Figure 2 :

 2 Figure 2: A C-arm with its flat-panel detector (dark rectangle, top left) and the Xray source on the opposite side of the C. The C-arm rotates around the patient bed. Left: The Patient Coordinate System (PCS). Middle: The Detector Coordinate System (DCS) and the corresponding 3D Detector Coordinate System (3DCS). w = u × v . Right: The Oriented Coordinate System (OCS) only differs from the 3DCS by the origin (O instead of O d ).

  Figure 5

Figure 6 :

 6 Figure 6: Top left: The central slice of a spherical object µ. Top right: The profile corresponding to the dotted red line. Middle: Backprojection of 3 (left) and 5 (right) projections. Bottom left: The central slice of the complete backprojection image (computed from 180 equally spaced projections). The sphere edges are blurred. Bottom right: The corresponding profile.

Figure 7 :

 7 Figure 7: The fanbeam geometry of one projection. The ray passing through x intersects the detector at u = r x • u/(rx • w ).

2 .

 2 Compute line integrals in the u-direction, 3. Differentiate the obtained 1D signal with respect to v, 4. Divide by cos 2 θ.

Figure 8 :

 8 Figure 8: Left: The 3D Detector coordinate system. Right: the spherical coordinates of the vector α. Note the unusual orientation: the polar axis is horizontal. The vector β has spherical coordinates ( β, 0).

Figure 9 :

 9 Figure 9: View of the plane u = 0.

Figure 10 :

 10 Figure 10: Given the PCS (O, e x , e y , e z ); the source s λ has coordinates (x λ , y λ , z λ ).In the 3DCS (O d , u, v , w ), the source has coordinates (u 0 , v 0 , f ).

Figure 11 :

 11 Figure 11: The zxy Tait-Bryan angles convention. First, a rotation of φ about the y world axis is applied (left). Then, a rotation of θ about the rotated x axis (middle). Finally, a rotation of η about the rotated z axis (right). The final coordinates are the OCS coordinates (x o , y o , z o ).

  x λ : x -coordinate of the source y λ : y -coordinate of the source z λ : z -coordinate of the source η : in-plane angle θ : out-of-plane angle φ : out-of-plane angle f : focal distance ( f > 0) u 0 : u -coordinate of principal point v 0 : v -coordinate of principal point

Figure 12 :

 12 Figure 12: The intrinsic matrix. The 3D point (x, y, z) is projected at the position (u, v) on the detector.Left: Top view (y = 0). Right: Side view (x = 0).

Figure 13 :

 13 Figure 13: The 2D Shepp-Logan phantom (left) and its sinogram over 180 degrees (right).

2 .

 2 Elements of S are denoted s and are represented by the triple r , θ, b where r > 0, θ ∈ [0, 2π[ and b ∈ R 2 . The non-negative real r is the scaling factor, θ is the rotation angle and b is the translation. The composition of two direct similarities is also a similarity. If R θ denote the rotation with angle θ, the transformation s = (r , θ, b) transforms x in s( x) with s( x) = r R θ x + b (114) The composition of s = r , θ, b (first) and s = r , θ , b (second), denoted ss , is easily found to be the similarity (r r , θ + θ , r R θ b + b ). Of course, the identity I belongs to S (the corresponding triple is (1, 0, 0)) and for s = (r , θ, b) ∈ S, the reciprocal transformation is

  Figure 15

Figure 17 :

 17 Figure 17: Fan-beam projections acquired along a line. The source is displaced on the vertical line x = 0. One ray is identified with the unit vector α or equivalently by the angle φ it makes with the horizontal axis (Ox).

Figure 19 :Figure 20 :

 1920 Figure 19: The integral of the weighted projections along lines parallel to the line between two source positions are equal.

Figure 21 :

 21 Figure 21: a. Left: The 3D Detector coordinate system. a. Right: the spherical coordinates of the vector α. Note the unusual orientation: the polar axis is horizontal. The vector β has spherical coordinates ( β, 0). b. View of the u = 0 plane.

  et al. "Two cone-beam consistency conditions for a circular trajectory." In: Proceedings of the fourfth international conference on image formation in x-ray computed tomography. Bamberg, Germany, 2016, pp. 431-434. We recall the intermediate function G τ introduced in Equation 55, Chapter I-1, for 0 ≤ τ ≤ 1,

Figure 22 :

 22 Figure 22: Circular cone-beam geometry. For each projection, the detector is aligned: the source, the origin of the PCS and the origin of the DCS are aligned and the vertical axis of the detector coincides with the axis of rotation y.

Figure 23 :

 23 Figure 23: Simulated truncated data. Three projections out of 72 were truncated as indicated. The dotted line is the intersection of the trajectory plane with the detector. Truncation of projection 10 (top right) is an axial truncation, away from the line. Truncation of projection 40 (bottom left) impacts the line. In projection 55 (bottom right), an interior truncation is simulated.

S

  I M U L AT I O N O F A N O N -R I G I D M O T I O N In the second experiment, one of the ellipsoids (a sphere of radius 4 mm) of the Shepp-Logan phantom was undergoing a vertical motion. The displacement started at projection index 20 and terminated at projection index 65. The trajectory of the displaced sphere intersected the central line. See Figure 24.

Figure 24 :

 24 Figure 24: Simulation of a non-rigid motion. A sphere of radius 4 mm is displaced, starting in projection 47, traversing the central line, and ending at projection 65.

Figure 25 :

 25 Figure 25: |Plots of the function Ḡ (top) and Ḡs (bottom) vs. the projection index in the case of truncation (left) and motion (right). The dotted line in the plot of the right column indicates the mean of Ḡ over the consistent projections. It helps appreciating the inconsistency (though it has no physical or theoretical significance).

Figure 26 :

 26 Figure 26: Same simulations as Figure 25 except that the line integral signals were smoothed with a smoothing kernel of size 11. See Equation 162.

Figure 27 :

 27 Figure 27: Evaluation of Ḡm (λ) for m = 0, ..., 5 for 50 projections, equally-spaced over a 360°scan. The black solid line corresponds to the trigonometric polynomial fitting of the data. The fitting problem was solved using a elementary least-squares approach. The outliers in projection indices 20 and 40 correspond to the inconsistency induced by the truncation of the data.

1 m = 2 Figure 28 :

 1228 Figure 28: The function Ḡm (λ) for m = 0, 1, 2 for 50 projections, equally-spaced over a 360°scan. The black solid line is a trigonometric polynomial, which bets fits the data in a least-squares sense.

  , v)ρ m (v)u m e i l λ du dv dλ = 0, |l | > m (184) or, by changing from the (u, v)-coordinates to the spherical coordinates in the 2D inner integral 2π 0 g λ α θ,φ cos φ ρ m (sin θ) tan m φ dφ dθ e ±i l λ dλ = 0, l > m ≥ 0. (185)

α

  = (sin ϕ cos θ, sin ϕ sin θ, cos ϕ) (186) with ϕ ∈ [0, 2π[ and θ ∈ [0, π[. The resulting consistency conditions assert that e ±2i nθ sin ϕ dϕ dθ e ±i l λ dλ = 0, (187)

Figure 29 :

 29 Figure 29: Picture of the CT system. Source (left) is fixed. Turntable (middle) and detector (right) are adjustable. Detector size: 35x35 mm. Pixel size: 17.09 µm.

Figure 30 :

 30 Figure 30: Illustration of the eight geometric parameters. The detector orientation is defined by three Euler angles. η is the in-plane angle. φ and θ are outof-plane rotations about the v-axis and the u-axis respectively.

Figure 32 :

 32 Figure 32: View of one plane P i , j ,k . Order-0 fanbeam DCCs state that the integral of the cosine-weighted projections are equal.

for

  Each pair of sources (i , j ) ∈ Ω: do 4:

6 :for 7 : 8 :

 678 Each row k: doCompute the line integrals G i , j ,k and G j ,i ,k .Compute the squared difference c i , j ,k .

Figure 34 :

 34 Figure 34: The cost function as a 1D function of each estimated parameter, evaluated on the simulated projection data of a standard Shepp-Logan phantom. 90 equally spaced projections were simulated over a full 360 • circular acquisition. For each parameter, the cost function is evaluated at 50 equally-spaced parameter values, ranging from -2 to 2. Note that the minimum function value is not zero due to numerical errors.

Figure 35 :

 35 Figure 35: Top and middle rows : Full-scan results. Top left: The 3D Shepp-Logan phantom. Top right: Reconstruction with true geometry. Middle left: Reconstruction with nominal geometry. Middle right: Reconstruction with DCC-calibrated geometry. Bottom left: Short scan reconstruction with nominal geometry. Bottom right: Short scan reconstruction with DCCcalibrated geometry.

Figure 36 :

 36 Figure 36: Consistency metric as a function of θ and v 0 (Top: 3D plot of the function. Bottom: gray-value 2D map and contour lines). All 3 others geometric parameters were set to their true values. The white solid line shows the valley along which inconsistency is minimal.

Figure 37 :

 37 Figure 37: Top row: Pictures of the imaged objects. Bottom row: One projection of each object.

Figure 38 :

 38 Figure 38: Coronal (left) and transverse (right) slices of the reconstructed image without calibration (top row), with DCC-based calibration (second row) and with marker-based calibration (third row). The intersection of both slices is represented by the white line. The corresponding intensity profiles are plotted on the bottom figure.

Figure 39 :

 39 Figure 39: Axial truncation management: only those rows between the two dashed lines are retained in the virtual projection.

Figure 40 :

 40 Figure 40: Concrete sample. Coronal (left) and transverse (right) slices of the reconstructed volume without calibration (top row) and with our DCCbased calibration (bottom row).

Figure 41 :

 41 Figure 41: Sponge sample. Coronal (left) and transverse (right) slices of the reconstructed volume without calibration (top row) and with our DCC-based calibration (bottom row).

Figure 42 :

 42 Figure 42: The cost function along the valley line for θ = 0.2 ± 10 degrees.

Figure 43 :

 43 Figure 43: Side view of the (y, z)-plane. The black thick lines represent the detectors. In its nominal position, the detector makes an angle θ with the rotation axis (O y). The centre of the detector O d , the origin O, and the source are on the same line. This corresponds to a detector offset v 0 .If the out-of-plane angle is modified to θ , the ray from s λ to the centre of the detector is displaced accordingly, introducing inconsistency. The inconsistency is counteracted by adjusting v 0 so that the ray returns to its original position.

Figure 44 :

 44 Figure 44: A pair of projections in general positions. The unmarked black arrows are the unit vectors of the PCS. The resulting epipolar geometry is characterized by the ECS (O, e e x , e e y , e e z ). e e x is parallel to the baseline b. The vector e e y is the normal to the plane containing O and the two source positions s 1 and s 2 .

Figure 45 :

 45 Figure 45: Spherical coordinates (θ, φ) are defined with respect to the epipolar geometry of the pair of projections. The baseline is the polar axis. The positions of the source in the ECS are (α 1 , 0, d ) and (α 2 , 0, d ) respectively. The point (u, v) lies on the virtual detector. We have the relations v = d tan θ and u = α 1 + d tan φ/ cos θ.

Figure 46 :

 46 Figure 46: The epipolar geometry of one pair of projections. All epipolar planes E θ intersect the detectors in epipolar lines l θ,i with i = 1, 2. All epipolar lines meet in the epipoles ε i .

Figure 47 :

 47 Figure 47: For an epipolar plane, the line v = v θ on the virtual plane (red) is projected on the physical detector to the epipolar line l θ .

S T E P 3 :

 3 H A N D L I N G T H E D I R A C δ D I S T R I B U T I O N . The term in G f with the Dirac function requires a particular treatment. It has the form δ (m(U ,V )) with

  , the remaining term | det H |/(d |∆ 13 |) is easily found to be 1/| cos( e o z , e e x )|.

3. 4 G

 4 R A N G E AT B A S E D C O N S I S T E N C Y C O N D I T I O N S A N D T H EA N G U -L A R F O R MWe remind theGrangeat intermediate function (see Equation 67 in page 20 and Figure 51 left)

Figure 51 :

 51 Figure 51: Two ways of computing G g ,i β . Left: The u-axis is parallel to the line of intersection. The weighted-projection is integrated in the u-direction, then differentiated in the v-direction. This is the original method, introduced by Grangeat. Right: The differentiation with respect to the angle of the epipolar line on the detector reflects the epipolar nature of a projection pair.

δ

  t α(γ) • β(θ) gλ e + t α(γ) t dt dγ (is the integral of the projection weighted by the cosine of the incidence angle and the inverse of the distance to the epipole. We let β D denote the orthogonal projection of β(θ) onto the detector plane:β D = β(θ) -β(θ) • w w . We have α(γ) • β(θ) = α(γ) • β D .Let γ θ denote the angle of the epipolar line l θ . The vector α(γ θ ) is perpendicular to β(θ) (because it lies on the epipolar plane E θ ) and to w (because it lies on the detector plane). We then have β D • α(γ θ ) = 0. Since α(γ θ ) = (cos γ θ , sin γ θ ), we have β D = ± β D (-sin γ θ , cos γ θ ) and (286)

2 = β D 2 + β(θ) • w 2 w = β D 2 + 2 = 53 2 =

 22222532 Figure 53

Figure 54 :

 54 Figure 54: The two projection images of the thorax Forbild phantom.

θFigure 55 :

 55 Figure 55: Top: The functions G f ,1 (θ) (green curve) and G f ,2 (θ) (blue curve) for θ ∈ [-36 • , 36 • ] (left). A zoom on the dashed box (right). Bottom: The same curves with G g ,i (θ) with i = 1, 2.

Figure 56 :

 56 Figure 56: Same plots as Figure 55 with noisy projections.

θFigure 57 :

 57 Figure 57: The plot on the right is a zoom on the dashed box of the plot on the left.Both curves coincide with very high accuracy.

  is compact. The set of test functions is denoted D(R) or simply D.The first non-trivial fact about D is that it is not restricted to the null function. The function φ defined by φ : x → e function. Indeed, φ is compactly supported on [-1, 1]. The derivatives of exp -1 1-x 2 are obtained by multiplying exp -1 1-x 2 by a rational function. As |x| → 1, such product tends to 0, hence the smoothness of φ.

) 1 . 1 . 2

 112 The space of distributions D By definition, a distribution is an element of the topological dual space of D i. e. the set of all continuous linear forms on D. This set is denoted D .

Definition 11 .

 11 Let T ∈ D a distribution. The derivative of T , denoted DT is, by definition< DT , φ >= -< T , φ > . (312)It is clear that DT is a linear form on D. Again, let φ k k -→ 0 and a compact K as in Definition 8. By definition of convergence in D, the convergence of φ k to zero induces the convergence of the sequence φ k to zero. Note that Suppφ

1. 1 A

 1 S H O R T I N T R O D U C T I O N T O D I S T R I B U T I O N S 123

) 1 . 1 . 5

 115 343) ρ(x) = |σ| e 2i πxσ dσ. (344Focus on the Dirac δ distribution

Rδ 3 . 24 (

 324 (m(x)) φ(x) dx = R R e 2i πσxτ(x) dσ φ(x) dx (of calculations leads to the following proposition, which is needed in Chapter II-Proposition Composition of δ with a smooth function). Let m be a smooth 1D function, with one simple root: m(x) = 0 if and only if x = 0 and m (0) = 0. Let φ a test function. Thenδ (m(x)) φ(x) dx = m (0)φ(0)m (0)φ (0) m (0) 2 |m (0)| . (352)Proof. We only sketch the proof.δ (m(x)) φ(x) dx = m (x)δ (m(x)) δ (m(x)) φ (x)m (x) -φ(x)m (x) 23, we finally obtain δ (m(x)) φ(x) dx = -φ (0)m (0) -φ(0)m (0) m (0) 2 |m (0)|(358)In the particular case of a function m satisfying m(0) = 0 and m (0) = 1, one gets δ (m(x)) φ(x) dx = -φ (0),

S 1 f α d α = 2π 0 fff

 10 (cos θ, sin θ) dθ, (360) where θ is the angle on the unit-circle. The proof of this proposition is trivial. Proposition 26 (Integration over S 2 ). With latitude-longitude spherical coordinates (θ, ϕ) ∈ [0, 2π[ × [-π/2, π/2[ of the sphere: cos θ cos ϕ, sin θ cos ϕ, sin ϕ cos ϕ dθ dϕ, (361) With colatitude-longitude coordinates (θ, ϕ) ∈ [0, 2π[ × [0, π[ of the sphere: cos θ sin ϕ, sin θ sin ϕ, cos ϕ sin ϕ dθ dϕ, (362)

Figure 58 :

 58 Figure58: A projective plane is the set of lines passing through the origin of a 3D vector space. By fixing a coordinate system, it is an affine plane plus points at infinity. The set of points at infinity (the horizon) are in a 1-1 correspondence with points on the semi-unit-circle in Π 0 .

Figure 59 :

 59 Figure 59: A representation of the projective plane. The bold line represents the points at infinity. All parallel lines (in the affine plane) meet at the same point at infinity.

Figure 60 :

 60 Figure 60: Geometry of one projection. The green rectangle on the left is the detector. The source position is s. All 3D points in the plane parallel to the detector and containing the source s are projected on points at infinity. On the contrary, any point which is not in this plane is projected on a finite point of the detector.

Figure 61 :

 61 Figure 61: The operation of resampling a projection into a virtual detector geometrically consists in a homography.

BProposition 28 .

 28 A P P E N D I X B : A U X I L IA R Y P R O O F S 2.1 H O M O G E N E O U S F U N C T I O N S O F D E G R E E -2In this appendix, we prove the following The set of 1D distributions which are homogeneous of degree -2 form a vector space of dimension 2, spanned by δ and ρ.Proof. Let H 2 = f ∈ D(R)s.t. ∀a > 0, x ∈ R, f (ax) = f (x)/a 2 and f ∈ H 2 .By taking the Fourier transform, we have, for a > 0, f (aσ) = f (s) e -2i π(aσ)s ds (376) = f (s /a) e -2i πσs ds/a (377)= a 2 f (s ) e -2i πσs ds/a (378)= a f (σ),

2. 4 C

 4 O M P U TAT I O N O F ∂γ/∂θ .We use the notation and the mathematical context of Section 3.5, where the Lemma 7 was introduced. We prove that

  definition, (0, cos θ, sin θ) T are the coordinates of β(θ) in the ECS. Left multiplying by R 12 gives the coordinates in the PCS. Left multiplying by R p gives the coordinates in the OCS and left multiplying by 1

  More precisely,

	Hounsfield built the
	first scanner in 1971,
	based on the
	theoretical works of
	Cormack, published
	in the 1950s.

  and placed at the origin O

	6 By aligned, we
	mean that 1. the
	principal ray passes
	through the centre of
	the trajectory O and
	the detector origin
	O d and 2. the
	direction of the
	detector is parallel
	to the x-axis when
	λ = 0.

Table 1 :

 1 Properties of the filters.

		Hom. deg.-1 Hom. deg. -2
	Even	δ	ρ
	Odd	h	δ

Table 2 :

 2 The nine geometric parameters of one projection.

  et al. "Data Consistency Conditions for Cone-Beam Projections on a Circular Trajectory." In: IEEE Signal Processing Letters 23.12 (2016), pp. 1746-1750.

  et al. "Calibration for Circular Cone-Beam CT Based on Consistency Conditions." In: IEEE Transactions on Radiation and Plasma Medical Sciences 1.6 (2017), pp. 517-526. ISSN: 2469-7311

Table 3 :

 3 Results on Simulated data

	η	θ	φ	u 0	v 0
	in deg. in deg. in deg. in mm. in mm.

(right). The set Ω was built with nine equally spaced projections over the 220 • angular range.

Table 4 :

 4 Calibration using different subsets Ω of projections

		η	θ	φ	u 0	v 0
		in deg. in deg. in deg. in mm. in mm.
	DCC-based #1 0.9301 0.2032 -0.1107 -0.0016 0.2912
	DCC-based #2	0.926	0.2198 -0.1091 -0.0013 0.2754
	DCC-based #3 0.9092 0.2439	-0.104 -0.0024 0.2846
	DCC-based #4 0.8844 0.2476 -0.1061 -0.0047 0.2997
	DCC-based #5 0.8567 0.2406 -0.1137 -0.0087 0.3108
	DCC-based #6 0.8542 0.2237 -0.1014 -0.0083	0.26
	DCC-based #7 0.8501 0.2075 -0.1099 -0.0025 0.2535
	DCC-based #8 0.8485 0.2159 -0.1097 -0.0005 0.3015
	SD	0.035	0.017	-0.004	-0.003	0.020

the edges of the object (see top-row of Figure

  [START_REF] Clackdoyle | Full data consistency conditions for cone-beam projections with sources on a plane[END_REF] = H 21 H 32 -H 22 H 31 . In order to factor out |∆ 13 | in the above equation, we shall prove that ∆ 13 = 0. By a well-known property of the inverse of a matrix, we have ∆ 13 = det H H-1 31 . The matrix H is easily inverted as

  n ≥ 2 : fp(1/x n ) We want to generalize the distribution pv (1/x) to function of the form 1/x n for n ≥ 2. The problem is that the nice behaviour with 1/x does not hold any more. Let f denote the function 1 R\[-, ] (x)/x 2 and T f the corresponding regular distribution. For all > 0 and |x| > one can write φ

Table 5 :

 5 Derivatives of usual distributions.

Table 6 :

 6 Fourier transform of usual distributions



  H 11 H 11 H 11 H 21 H 22 H 23 H 31 H 32 H 33

  H 11 H T 3• X -H 31 H T 1• X H 22 H T 3• X -H 32 H T 2• X -(H 11 H 22 -H 21 H 12 ) + H T 1• X H T 3• X (-H 31 H 22 + H 21 H 32 ) + X (H 11 H 22 -H 21 H 12 ) + H T 2• X (-H 11 H 32 + H 31 H 12 ) + H 11 H 12 H 11 U + H 12 V + H 13 H 21 H 22 H 21 U + H 22 V + H 23 H 31 H 32 H 31 U + H 32 V + H 33

	The Jacobian determinant is Jac = ∂u ∂U	∂v ∂V -∂u ∂V	∂v ∂U . Easy calculations give
	Jac =	1 H T 3• X
			H 21 H T 3• X -H 31 H T 2• X H 12 H T 3• X -H 32 H T 1• X	(385)
	=	1 H T 3• X	4 H T 3• X
					H T 2• X H T 3• X (-H 11 H 32 + H 31 H 12 )	(386)
	=	1 H T 3• X	3 H T 3• H T 1• X (-H 31 H 22 + H 21 H 32 )	(387)
	=	1 H T 3• X	3 det	   	H 11 H 12 H T 1• X H 31 H 32 H T 3• X H 21 H 22 H T 2• X	   	(388)
								
	=	1 H T 3• X	3 det	  				  
								133

4 2

  The angle γ θ of the epipolar line l θ is given by Since cos 2 γ θ = b 2 /(a 2 + b 2 ), one finally gets

	tan γ θ =	-a b	=	-H 21 + v θ H 31 H 22 -v θ H 32	.	(404)
	Straight forward computations lead to
	∂ ∂θ	tan γ θ =	d cos 2 θ	-∆ 13 (H 22 -v θ H 32 ) 2 .	(405)
	Using Equation 267, we obtain
	∂ ∂θ	tan γ θ =	d 2 f 2 cos 2 θ	cos( e o z , e e x ) (H 22 -v θ H 32 ) 2 =	cos( e o z , e e x ) cos 2 θ	a 2 + b 2 b 2 .	(406)
	Finally,
	∂ ∂θ	tan γ θ =	∂γ θ ∂θ	×	∂ ∂γ	tan γ θ	(407)
				=	∂γ θ ∂θ	×	1 cos 2 γ θ	.	(408)
	∂γ ∂θ	=	cos( e o z , e e x ) cos 2 θ	,	(409)
	which ends the proof.

T O M O G R A P H Y : A N O V E R V I E W

Note that the

See https://en.wikipedia.org/wiki/Euler _ angles for detailed presentation of Euler angles.