# Cone-beam pair-wise data consistency conditions in helical CT

Mélanie Mouchet, Jean Michel Létang, Jérome Lesaint, Simon Rit

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mélanie Mouchet, Jean Michel Létang, Jérome Lesaint, Simon Rit. Cone-beam pair-wise data consistency conditions in helical CT. IEEE Transactions on Medical Imaging, 2023, 42 (10), pp.2853-2864. 10.1109/TMI.2023.3265812. hal-04072471

HAL Id: hal-04072471

https://hal.science/hal-04072471

Submitted on 18 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cone-beam pair-wise data consistency conditions in helical CT

Mélanie Mouchet, Jean Michel Létang, Jérôme Lesaint and Simon Rit

**Abstract**— Data consistency conditions (DCC) are mathematical equations characterizing the redundancy in X-ray projections. They have been used to correct inconsistent projections before computed tomography (CT) reconstruction. This article investigates DCC for a helical acquisition with a cylindrical detector, the geometry of most diagnostic CT scanners. The acquired projections are analyzed pair-by-pair. The intersection of each plane containing the two source positions with the corresponding cone-beams defines two fan-beams for which a DCC can be computed. Instead of rebinning the two fan-beam projections to a conventional detector, we directly derive the DCC in detector coordinates. If the line defined by two source positions intersects the field-of-view (FOV), the DCC presents a singularity which is accounted for in our numerical implementation to increase the number of DCC compared to previous approaches which excluded these pairs of source positions. Axial truncation of the projections is addressed by identifying for which set of planes containing the two source positions the fan-beams are not truncated. The ability of these DCC to detect breathing motion has been evaluated on simulated and real projections. Our results indicate that the DCC can detect motion if the baseline and the FOV do not intersect. If they do, the inconsistency due to motion is dominated by discretization errors and noise. We therefore propose to normalize the inconsistency by the noise to obtain a noise-aware metric which is mostly sensitive to inconsistencies due to motion. Combined with a moving average to reduce noise, the derived DCC can detect breathing motion.

**Index Terms**— data consistency conditions, helical CT, motion detection, statistical noise.

## I. INTRODUCTION

DATA consistency conditions (DCC) are mathematical equations which characterize the range of a linear operator, i.e. in computed tomography (CT), the range of the Radon transform for the 2D parallel geometry or the

This work was supported in part by the French National Association for Technical Research (ANRT), in part by SIRIC LYriCAN under Grant INCa-INSERM-DGOS-12563, and in part by the Labex PRIMES of Université de Lyon, within the Program Investissements d'Avenir operated by the ANR under Grant ANR-11-LABX-0063 and ANR-11-IDEX-0007. (Corresponding author: Mélanie Mouchet.)

The patient data included in this study comply with current regulations for human protection and GDPR regulations (MR004).

Mélanie Mouchet is with Siemens Healthcare, Paris, France and Univ Lyon, INSA-Lyon, Université Claude Bernard Lyon 1, UJM-Saint Étienne, CNRS, Inserm, CREATIS UMR 5220, U1294, F-69373, Lyon, France (e-mail: melanie.mouchet@yahoo.com).

Jean Michel Létang, Jérôme Lesaint, and Simon Rit are with Univ Lyon, INSA-Lyon, Université Claude Bernard Lyon 1, UJM-Saint Étienne, CNRS, Inserm, CREATIS UMR 5220, U1294, F-69373, Lyon, France (e-mails: lesaint.jerome@gmail.com, jean.letang@insa-lyon.fr, simon.rit@creatis.insa-lyon.fr).

divergent beam transform for divergent geometries [1]. In other words, DCC characterize the redundancy between x-ray projections and must be verified by the projections. For example, the parallel projections of a 2D function should verify the following consistency condition: the sum of the pixel values of each projection is the same for all measured projections. If it is not verified, the projections cannot be the Radon transform of some attenuation function and the projections are *inconsistent*. This DCC is known as the zeroth order of the Helgason-Ludwig consistency conditions [2], [3]. Practical applications in CT are numerous because DCC provide mathematical relationships to correct the projections acquired by the CT scanner, prior to tomographic reconstruction. The correction of artifacts in CT using the DCC is based on the minimization of a cost function derived from the DCC. Past examples include the completion of missing projections in limited angle tomography [4], the correction of ring artifacts due to malfunctioning detectors [5], the correction of beam hardening in polychromatic CT [6]–[9], the calibration of geometric parameters such as the detector position and orientation [10]–[12], or the detection and the compensation of rigid patient motion [13]–[17].

DCC are specific to the geometry of the CT scanner since the geometry defines the lines along which the linear operator integrates to model the x-ray acquisition process. For parallel projections of a 2D function, complete (i.e. necessary and sufficient) DCC are known as Helgason-Ludwig conditions [2], [3]. For fan-beam projections, various DCC have also been proposed, e.g., complete DCC for sources on a line [18]. The latter have been derived from the Helgason-Ludwig DCC and have a similar expression. Those DCC state that the moments of order *n* of the parallel or fan-beam projections are a homogeneous polynomial of order *n* in the trajectory coordinate. In this work, our focus is on cone-beam projections for which some DCC have been developed for specific source trajectories. The cone-beam equivalent of the Helgason-Ludwig DCC for sources on a trajectory plane which does not intersect the object was derived in [19]. Necessary DCC for sources on a circular trajectory where the trajectory plane intersects the object are described in [20]. These DCC use an intermediate function which should be constant over the entire trajectory. By studying the link between [19] and [20], other DCC for other geometrical configurations of trajectory and detector planes have been derived in [21].

Other DCC can be used for any source trajectory by
considering pairs of cone-beam projections. Those DCC state that some intermediate function of the projection values should be equal for the two projections of the pair. Such DCC have been derived from Grangeat's relation [22] and provide a set of DCC for each source pair, one for each plane containing the line defined by the corresponding source positions [10], [11], which we refer to as the epipolar plane and the baseline as in [11]. Each projection is first weighted by the cosine of the incidence angle and the derivative is taken in the direction perpendicular to the epipolar line, i.e. the line of intersection between the projection and the epipolar plane. DCC state that the sums of these pre-processed projections along corresponding epipolar lines are equal. A similar pair-wise approach was followed by Lesaint et al [12] considering two fan-beam projections for each epipolar plane and applying fan-beam DCC [15], [23], which are the zeroth order DCC of more complete DCC defined in [18]. Those zeroth order DCC, also referred to as integral invariants, have been known for quite some time and derived in many ways: in the Fourier domain in [5], from some homogeneous function and its symmetric group in [13] or from John's equation in [15]. In [24], they were all proven mathematically equivalent and it was later proven that if these pair-wise fan-beam DCC are verified, then the Grangeat-based pair-wise DCC will also be verified [25].

Diagnostic CT scanners generally have an x-ray source that follows a helical trajectory coupled with a cylindrical detector. This article addresses DCC for this geometry for which, to our knowledge, DCC have not been investigated so far. A direct approach is to use pair-wise DCC by rebinning the projections acquired by a cylindrical detector to projections on a plane parallel to the baseline, as proposed in [12]. Such an approach would require that the baseline does not intersect the object and we propose instead to compute them on the cylindrical detectors directly to alleviate this condition and increase the number of DCC. This induces a singularity in the numerical implementation of the DCC which we address. Another characteristic of diagnostic scanners is the limited size of the detector in the axial direction which entails axial truncation of the cone-beam projections. Two cone-beam projections may scan completely different parts of the object and DCC are only available for a subset of pairs of projections. This article provides the conditions that two source positions must verify to be part of this subset if the projections are only truncated axially. We also derive the DCC uncertainty stemming from the noise of the projections to implement a noise-aware consistency metric. The metric is validated using numerical simulations and real data of a phantom and a patient for the detection of breathing motion.

## II. MATERIALS AND METHODS

### A. Description of the geometry

The geometry of the helical scanner is illustrated in Fig. 1. The three-dimensional (3D) coordinates $(x, y, z)$ of the scanned object or patient are defined in the orthogonal basis $(e_x, e_y, e_z)$. The source moves along the helix

<figure><img src="image_3.png" /><figcaption>Fig. 1: Cone-beam geometry of the helical scanner used in this study.</figcaption></figure>

$s_λ = (R \cos λ, R \sin λ, dλ/2π)$, where $λ$ is the gantry angle, taken between $e_x$ and the projection of $s_λ$ on the $(e_x, e_y)$ plane (measured positively in the clockwise direction when viewed in the direction $e_z$), and $d$ is the helix pitch, i.e., the distance traveled by the source along the z-axis in one rotation.

The two-dimensional (2D) detector is a portion of cylinder of radius $D$ which axis is defined by the source point $s_λ$ and the rotation axis $e_z$. At a given source position $s_λ$, the 3D coordinates of a detector point with coordinates $(γ, v)$ are

$$ \mathbf{p}_{\lambda}(\gamma, v) = s_{\lambda} + \begin{pmatrix} -D \cos(\gamma - \lambda) \\ D \sin(\gamma - \lambda) \\ v \end{pmatrix}. \quad (1) $$

The angle $γ$ between the x-ray lines at $p_λ(0,0)$ and $p_λ(γ,0)$ is measured positively in the counterclockwise direction when viewed in the direction $e_z$. The detector limits are $γ \in [-γ_{max}, γ_{max}]$ (with $γ_{max} < π/2$) and $v \in [-v_{max}, v_{max}]$. The field-of-view (FOV) is the region defined by $\sqrt{x^2 + y^2} \le R_{FOV}$, with $R_{FOV} = R \sin γ_{max}$.

The scanner probes the linear attenuation coefficient $f$ of the object in the FOV according to the line integral model

$$ g_{\lambda}(\gamma, v) = \int_{0}^{\infty} f \left( s_{\lambda} + t \frac{p_{\lambda}(\gamma, v) - s_{\lambda}}{\sqrt{D^2 + v^2}} \right) dt. \quad (2) $$

We assume that $f$ is zero outside the FOV, i.e., that projections $g_λ$ are not truncated laterally (but can be axially).

### B. Pair-wise geometry

The DCC used in this work are derived from fan-beam DCC for sources on a line as in [15], [23]. Any pair of source positions $(s_λ, s_{λ'})$ defines such a baseline $B_{λ,λ'}$ oriented with the unit vector

$$ \mathbf{b}_{\lambda,\lambda'} = \mathrm{sign}(\lambda' - \lambda) \frac{s_{\lambda'} - s_{\lambda}}{\|s_{\lambda'} - s_{\lambda}\|} \quad (3) $$

illustrated in Fig. 2.
<figure><img src="image_4.png" /><figcaption>Fig. 2: 3D illustration of the cone-beam pair-wise geometry. The baseline $B_{\lambda,\lambda'}$ is defined by the source positions $s_{\lambda}$ and $s_{\lambda'}$ (with $\lambda < \lambda'$ in this case). The intersection of the reference epipolar plane $\Pi_{\lambda,\lambda'}$ with the two detectors defines the two curves $C_{\lambda,\lambda'}$ and $C_{\lambda',\lambda}$ (thick dotted lines). The intersection of an epipolar plane $\Pi_{\tau}$ with the two detectors defines the two curves $C_{\tau}$ and $C_{\tau'}$ with $\tau = (\lambda, \lambda', \beta)$ and $\tau' = (\lambda', \lambda, \beta)$ (thick continuous lines). The dotted limits of the detectors are below the reference epipolar plane $\Pi_{\lambda,\lambda'}$ and the dashed ones are below the epipolar plane $\Pi_{\tau}$.</figcaption></figure>

Fan-beam DCC may be computed in any epipolar plane, i.e., any plane which contains the baseline $B_{\lambda,\lambda'}$. The reference epipolar plane $\Pi_{\lambda,\lambda'}$ is defined by its orthogonal

$$n_{\lambda,\lambda'} = c_{\lambda,\lambda'} \times b_{\lambda,\lambda'} \quad (4)$$

where

$$c_{\lambda,\lambda'} = \cos \bar{\lambda} e_x + \sin \bar{\lambda} e_y \quad (5)$$

is the unit vector orthogonal to the baseline $B_{\lambda,\lambda'}$ and the rotation axis $e_z$. The vector $c_{\lambda,\lambda'}$ points from the midpoint between $s_{\lambda}$ and $s_{\lambda'}$ (which is on the baseline $B_{\lambda,\lambda'}$) towards the helix axis, with

$$\bar{\lambda} = \begin{cases} \frac{\lambda + \lambda'}{2} & \text{if } \frac{\lambda' - \lambda}{2} (\text{mod } 2\pi) \in [\frac{\pi}{2}, \frac{3\pi}{2}], \\ \frac{\lambda + \lambda'}{2} + \pi & \text{otherwise.} \end{cases} \quad (6)$$

The other epipolar planes $\Pi_{\tau}$ with $\tau = (\lambda, \lambda', \beta)$ are defined by their angle $\beta$ with the reference plane $\Pi_{\lambda,\lambda'}$ (Fig. 2). Their orthogonal is

$$n_{\tau} = \cos \beta n_{\lambda,\lambda'} - \sin \beta c_{\lambda,\lambda'}. \quad (7)$$

The angle $\beta$ is measured positively in the counterclockwise direction when viewed in the direction $b_{\lambda,\lambda'}$. With the third vector $c_{\tau} = b_{\lambda,\lambda'} \times n_{\tau}$, the triplet $(c_{\tau}, b_{\lambda,\lambda'}, n_{\tau})$ is an ortho-normal basis which is used in the following to derive cone-beam DCC from fan-beam DCC. If we let $c_{\tau} = (c_1, c_2, c_3)$, $b_{\lambda,\lambda'} = (b_1, b_2, b_3)$ and $n_{\tau} = (n_1, n_2, n_3)$ be the coordinates of these vectors, they

can be derived by three successive rotations

$$\begin{pmatrix} c_1 & b_1 & n_1 \\ c_2 & b_2 & n_2 \\ c_3 & b_3 & n_3 \end{pmatrix} = \begin{pmatrix} \cos \bar{\lambda} & -\sin \bar{\lambda} & 0 \\ \sin \bar{\lambda} & \cos \bar{\lambda} & 0 \\ 0 & 0 & 1 \end{pmatrix} \times \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \alpha & -\sin \alpha \\ 0 & \sin \alpha & \cos \alpha \end{pmatrix} \times \begin{pmatrix} \cos \beta & 0 & -\sin \beta \\ 0 & 1 & 0 \\ \sin \beta & 0 & \cos \beta \end{pmatrix} \quad (8)$$

with $\alpha = \arccos(b_{\lambda,\lambda'} \cdot (e_z \times c_{\lambda,\lambda'}))$ the angle between the axial planes (defined by $z = \text{const}$) and $B_{\lambda,\lambda'}$. We obtain

$$\left\{ \begin{aligned} c_1 &= \cos \bar{\lambda} \cos \beta + \sin \bar{\lambda} \sin \alpha \sin \beta, \\ c_2 &= \sin \bar{\lambda} \cos \beta - \cos \bar{\lambda} \sin \alpha \sin \beta, \\ c_3 &= \cos \alpha \sin \beta, \\ b_1 &= -\sin \bar{\lambda} \cos \alpha, \\ b_2 &= \cos \bar{\lambda} \cos \alpha, \\ b_3 &= \sin \alpha, \\ n_1 &= -\cos \bar{\lambda} \sin \beta + \sin \bar{\lambda} \sin \alpha \cos \beta, \\ n_2 &= -\sin \bar{\lambda} \sin \beta - \cos \bar{\lambda} \sin \alpha \cos \beta, \\ n_3 &= \cos \alpha \cos \beta. \end{aligned} \right. \quad (9)$$

### C. Fan-beam DCC in an epipolar plane

The intersection of the epipolar plane $\Pi_{\tau}$ and the cone-beam at source position $s_{\lambda}$ defines a fan-beam which may be parameterized by the angle $\phi$ between $c_{\tau}$ and each ray (Fig. 3). We note

$$g_{\tau}(\phi) = \int_{0}^{\infty} f(s_{\lambda} + t(\cos \phi c_{\tau} + \sin \phi b_{\lambda,\lambda'})) dt \quad (10)$$

the line integral at angle $\phi$ in the fan-beam. The zeroth order DCC for fan-beam projections along a line [18] states that the intermediate function, which is referred to a moment in the following,

$$M_{\tau} = \int_{-\pi}^{\pi} \frac{g_{\tau}(\phi)}{\cos \phi} d\phi \quad (11)$$

is constant for all source positions along the line. This integral, which presents two singularities at $-\pi/2$ and $\pi/2$, must be understood in the sense of the Cauchy principal value.

The computation of the moment $M_{\tau}$ requires untruncated fan-beam projections. Since the cone-beam projections are generally truncated in the $v$ direction, which is known as the long-object problem, we exclude projection pairs which define a baseline in the direction $e_z$ and epipolar planes with a normal $n_{\tau}$ orthogonal to $e_z$ by setting $\lambda' - \lambda \neq 0 (\text{mod } 2\pi)$, $\alpha \in [0, \pi) \setminus \{\pi/2\}$ and $\beta \in (-\pi/2, \pi/2)$. Consequently, the two source positions $s_{\lambda}$ and $s_{\lambda'}$ are the only ones along $B_{\lambda,\lambda'}$ and (11) only provides the following DCC in the epipolar plane $\Pi_{\tau}$

$$M_{\tau} = M_{\tau'} \quad (12)$$

where $M_{\tau'}$ is the moment of the projection $g_{\tau'}$ with $\tau' = (\lambda', \lambda, \beta)$ and $\beta$ in a range derived in (II-E).

We derive the expression of the moment $M_{\tau}$ in the detector coordinates $(\gamma, v)$. At source position $s_{\lambda}$, the
<figure><img src="image_5.png" /><figcaption>Fig. 3: Geometric illustration in the epipolar plane $\Pi_{\tau}$ of the angle $\phi$ between $c_{\tau}$ and the line defined by $s_{\lambda}$ and $p_{\lambda}(\gamma, v_{\tau})$. $\Pi_{\tau}$ is not orthogonal to the helix axis $e_z$</figcaption></figure>

epipolar plane intersects the detector along the curve $C_{\tau}$ defined by $(p_{\lambda}(\gamma, v_{\tau}) - s_{\lambda}) \cdot n_{\tau} = 0$, i.e.,

$$C_{\tau}: v_{\tau}(\gamma) = \frac{D}{n_3}(n_1 \cos(\gamma - \lambda) - n_2 \sin(\gamma - \lambda)) \quad (13)$$

where $n_3 = \cos \alpha \cos \beta \neq 0$ since $\alpha \in [0, \pi) \setminus \{\pi/2\}$ and $\beta \in (-\pi/2, \pi/2)$. The ray direction $p_{\lambda}(\gamma, v_{\tau}(\gamma)) - s_{\lambda}$ now depends on $\gamma$ only. We have a one-to-one mapping between the angles $\phi$ and $\gamma$ and we can change variable from $\phi$ to $\gamma$ in the integral defined by (11). The change of variable is available in the Appendix and we obtain the following expression for the projection moment in the cylindrical coordinates

$$M_{\tau} = \pi \int_{-\gamma_{\max}}^{\gamma_{\max}} \frac{\sin(\gamma_{\lambda, \lambda'}^*) h(\gamma_{\lambda, \lambda'}^* - \gamma)}{|\cos \alpha| \sin(\gamma_{\lambda, \lambda'}^* - \gamma)} \tilde{g}_{\lambda}(\gamma, v_{\tau}(\gamma)) d\gamma, \quad (14)$$

where $\tilde{g}_{\lambda}(\gamma, v_{\tau}(\gamma)) = Dg_{\lambda}(\gamma, v_{\tau}(\gamma))/\sqrt{D^2 + v_{\tau}^2(\gamma)}$, $h(\gamma) = 1/(\pi\gamma)$ is the Hilbert kernel and $\gamma_{\lambda, \lambda'}^*$ is a singularity corresponding to the angular coordinate at which the baseline may intersect the detector.

## D. Discretization

The detector has $N_{\text{cols}} \times N_{\text{rows}}$ pixels of size $\Delta_{\gamma} \times \Delta_{v}$. Their coordinates $\gamma$ and $v$ are $\gamma_i = -\gamma_{\max} + (i-1)\Delta_{\gamma}$, $i \in \{1, ..., N_{\text{cols}}\}$ and $v_j = -v_{\max} + (j-1)\Delta_{v}$, $j \in \{1, ..., N_{\text{rows}}\}$. The integral in (14) is discretized as follows:

$$M_{\tau} \approx \pi \sum_{i=1}^{N_{\text{cols}}} \frac{\sin(\gamma_{\lambda, \lambda'}^*) h_{\nu}(\gamma_{\lambda, \lambda'}^* - \gamma)}{|\cos \alpha| \sin(\gamma_{\lambda, \lambda'}^* - \gamma)} \tilde{g}_{\lambda}(\gamma_i, v_{\tau}(\gamma_i)) \Delta_{\gamma} \quad (15)$$

where $v_{\tau}(\gamma_i) \in [v_j, v_{j+1})$ is the v-coordinate of the point on the curve $C_{\tau}$ at angle $\gamma = \gamma_i$ such that $j = \lfloor (v_{\tau}(\gamma_i) + v_{\max})/\Delta_{v} \rfloor$ depends on $i$. The projection value is obtained with linear interpolation

$$g_{\lambda}(\gamma_i, v_{\tau}(\gamma_i)) = (1 - w_{i,j}) g_{\lambda}(\gamma_i, v_j) + w_{i,j} g_{\lambda}(\gamma_i, v_{j+1}) \quad (16)$$

<figure><img src="image_6.png" /><figcaption>Fig. 4: Example of a plane $\Pi_{\tau}$ for which the DCC cannot be computed because the curves $C_{\tau}$ and $C_{\tau'}$ do not verify the condition (18). The truncated part of the curves are shown in red.</figcaption></figure>

with $w_{i,j} = (v_{\tau}(\gamma_i) - v_j)/\Delta_{v}$ the linear interpolation weight. The kernel $h_{\nu}(\gamma)$ is a band-limited version of $h(\gamma)$. It is obtained by taking the inverse Fourier transform of the product between the Fourier transform of the Hilbert kernel and a Hann window. The resulting expression is

$$h_{\nu}(\gamma) = \frac{1 - \cos\left(\nu \frac{\pi\gamma}{\Delta_{\gamma}}\right)}{2\pi\gamma} + \frac{\gamma \left(1 + \cos\left(\nu \frac{\pi\gamma}{\Delta_{\gamma}}\right)\right)}{2\pi \left(\gamma^2 - \frac{\Delta_{\gamma}^2}{\nu^2}\right)} \quad (17)$$

where $\nu \in (0, 1]$ is a fraction of Nyquist's frequency and controls the bandwidth of the kernel.

## E. Applicability intervals for $\alpha$ and $\beta$

The DCC are only valid if the fan-beam projections $g_{\tau}(\phi)$ (10) are not truncated. Since the projections may be truncated in the axial direction, DCC are only applicable for some pairs of source positions (i.e. some $\lambda - \lambda'$) and some epipolar planes (i.e. some $\beta$) for which the detector curves $C_{\tau}$ and $C_{\tau'}$ are between the inferior and superior limits of the detector (see Fig. 2):

$$|v_{\tau}(\gamma)| \leq v_{\max}. \quad (18)$$

If this inequality is not verified for one of the two curves $C_{\tau}$ or $C_{\tau'}$, as illustrated in Fig. 4, the resulting fan-beam projections $g_{\tau}$ and $g_{\tau'}$ are truncated and the DCC cannot be used.

Replacing $(n_1, n_2, n_3)$ in $v_{\tau}(\gamma)$ (13) by their expression in (9) and $\lambda - \lambda'$ by $\gamma_{\lambda, \lambda'}^*$ from (41), we have, similarly to (43), $\sin(\gamma + \bar{\lambda} - \lambda) = -\sin(\gamma_{\lambda, \lambda'}^* \cos(\gamma - \gamma_{\lambda, \lambda'}^*))$ and the previous inequality becomes

$$D \left| \tan \alpha \cos(\gamma - \gamma_{\lambda, \lambda'}^*) + \frac{\tan \beta}{\cos \alpha} \sin(\gamma - \gamma_{\lambda, \lambda'}^*) \right| \leq v_{\max} \quad (19)$$

which gives the following conditions on the angles $\alpha$ and $\beta$:

* If $\gamma = \gamma_{\lambda, \lambda'}^*$ (which implies $\gamma_{\lambda, \lambda'}^* \in [-\gamma_{\max}, \gamma_{\max}]$), we have

$$|\tan \alpha| \leq \frac{v_{\max}}{D}. \quad (20)$$

This constraint on $\alpha$ indirectly constrains the interval separating two source positions where DCC can be applied
because $\alpha$ increases with the angular distance $|\lambda' - \lambda|$ between the pair of source positions. This specific condition where $\gamma$ coincides with the singularity, i.e. the direction of the baseline, does not depend on $\beta$ since $\beta$ is the angle defining the epipolar planes around the baseline.

Using $|\tan \alpha| = b_3/\sqrt{b_1^2 + b_2^2}$ yields

$$ \left| \frac{\lambda' - \lambda}{\sin\left(\frac{\lambda' - \lambda}{2}\right)} \right| \le \frac{4\pi R v_{\max}}{d D} \quad (21) $$

which shows that the smaller the helical pitch $d$, the larger the length of the applicability interval. The sine function in this inequality changes the applicability interval from rotation to rotation. Almost full rotations are applicable ((21) verified) for the rotations around the reference position $\lambda$ (i.e. when $|\lambda' - \lambda|$ is small with respect to the right-hand side of (21)). Only a part of the rotation will be applicable when hitting the limit, in intervals centered around values $(\lambda' - \lambda)$ mod $2\pi = \pi$ with a length decreasing when the angular distance to the reference $|\lambda' - \lambda|$ increases.

• If $\gamma \neq \gamma_{\lambda, \lambda'}^*$, then isolating $\beta$ from (19) yields

$$ l_{\lambda, \lambda'}(\gamma) \le \beta \le u_{\lambda, \lambda'}(\gamma) \quad (22) $$

where

$$ l_{\lambda, \lambda'}(\gamma) = \arctan \left( - \frac{\sin \alpha}{\tan(\gamma - \gamma_{\lambda, \lambda'}^*)} - \frac{v_{\max}}{D} \left| \frac{\cos \alpha}{\sin(\gamma - \gamma_{\lambda, \lambda'}^*)} \right| \right) \quad (23) $$

and

$$ u_{\lambda, \lambda'}(\gamma) = \arctan \left( - \frac{\sin \alpha}{\tan(\gamma - \gamma_{\lambda, \lambda'}^*)} + \frac{v_{\max}}{D} \left| \frac{\cos \alpha}{\sin(\gamma - \gamma_{\lambda, \lambda'}^*)} \right| \right) \quad (24) $$

are the equations of the smallest and highest angle $\beta$ for which the epipolar curve $v_{\tau}(\gamma)$ is not truncated for each detector column $\gamma$. Because $\gamma_{\lambda, \lambda'}^*$ is not the same for $s_{\lambda}$ and $s_{\lambda'}$, (22) must be verified for the two source positions. However, it can be shown that $l_{\lambda', \lambda}(\gamma) = -u_{\lambda, \lambda'}(-\gamma)$ and $u_{\lambda', \lambda}(\gamma) = -l_{\lambda, \lambda'}(-\gamma)$ because a $\pi$ rotation around the axis defined by $(s_{\lambda} + s_{\lambda'})/2$ and $c_{\lambda, \lambda'}$ swaps the two detectors. We therefore obtain the condition:

$$ |\beta| \le \beta_{\max} = \min_{\gamma \in [-\gamma_{\max}, \gamma_{\max}]} \min(u_{\lambda, \lambda'}(\gamma), -l_{\lambda, \lambda'}(-\gamma)). \quad (25) $$

If $\beta_{\max} \le 0$, the detector is too small to compute a DCC for this pair of source positions. Otherwise, we sample $\beta$ in $B = \lfloor 2\beta_{\max}D/\Delta_v \rfloor$ planes.

## F. Consistency metric

We define the consistency metric between two cone-beam projections as

$$ E_{\lambda, \lambda'} = \frac{B^{-1} \sum_{b=1}^{B} \left| M_{\tau_b} - M_{\tau_b'} \right|}{\sqrt{B \left( \text{Var}(\overline{M_{\lambda, \lambda'}}) + \text{Var}(\overline{M_{\lambda', \lambda}}) \right)}}. \quad (26) $$

The numerator of $E_{\lambda, \lambda'}$ is the mean absolute deviation from the DCC defined in (12) over the $B$ epipolar planes, which expected value is independent of $B$ if the samples are uncorrelated. The denominator is an approximate estimation of the noise (standard deviation) of the numerator without the absolute values, which assumes that $\overline{M_{\lambda, \lambda'}}$ and $\overline{M_{\lambda', \lambda}}$ are uncorrelated (to split the variance of the difference in a sum of variances) and uses the sum of the means of the $B$ moments

$$ \overline{M_{\lambda, \lambda'}} = \frac{1}{B} \sum_{b=1}^{B} M_{\tau_b}. \quad (27) $$

The $\sqrt{B}$ weight in the denominator of $E_{\lambda, \lambda'}$ makes it independent of $B$, i.e., compensates the decrease of $\text{Var}(\overline{M_{\lambda, \lambda'}})$ with the increase of $B$. The normalization by the noise (denominator) is necessary because the number of photons (and the resulting noise) strongly varies from projection to projection, e.g. between lateral and antero-posterior views of a patient. The impact of this noise on $M_{\tau}$ depends on the weight of each projection value $g_{\lambda}(\gamma, v)$ which also strongly varies with $\tau$ and $\gamma$ in the computation of the moments $M_{\tau}$ in (15).

The variance $\text{Var}(\overline{M_{\lambda, \lambda'}})$ can be derived following the approach used for fan-beam CT reconstruction in [26]. We assume that the noise between pixels of the projections is uncorrelated. The moments computed on different projections are uncorrelated, as previously assumed, but the $B$ moments computed on the same projection may use the same pixels and be correlated. Their co-variances should be accounted for with

$$ \text{Var}(\overline{M_{\lambda, \lambda'}}) = \frac{1}{B^2} \sum_{b=1}^{B} \text{Var}(M_{\tau_b}) + \frac{2}{B^2} \sum_{1 \le k < l \le B} \text{Cov}(M_{\tau_k}, M_{\tau_l}). \quad (28) $$

The terms in the sums above are expanded using (15). We note

$$ a_{\tau}(\gamma) = \Delta_{\gamma} \frac{\pi D}{\sqrt{D^2 + v_{\tau}^2(\gamma)}} \frac{\text{sign} \gamma_{\lambda, \lambda'}^* h_{\nu}(\gamma_{\lambda, \lambda'}^* - \gamma)}{|\cos \alpha| \text{sinc}(\gamma_{\lambda, \lambda'}^* - \gamma)} \quad (29) $$

and obtain

$$ \text{Var}(M_{\tau_b}) = \sum_{i=1}^{N_{\text{cols}}} a_{\tau_b}^2(\gamma_i)(1 - w_{i,j})^2 \text{Var}(g_{\lambda}(\gamma_i, v_j)) + \sum_{i=1}^{N_{\text{cols}}} a_{\tau_b}^2(\gamma_i)w_{i,j}^2 \text{Var}(g_{\lambda}(\gamma_i, v_{j+1})) \quad (30) $$

$$ \text{Cov}(M_{\tau_k}, M_{\tau_l}) = \sum_{i=1}^{N_{\text{cols}}} a_{\tau_k}(\gamma_i) a_{\tau_l}(\gamma_i) V_{i,j_k,j_l} \text{Var}(g_{\lambda}(\gamma_i, v_{j_k})) + \sum_{i=1}^{N_{\text{cols}}} a_{\tau_k}(\gamma_i) a_{\tau_l}(\gamma_i) W_{i,j_k,j_l} \text{Var}(g_{\lambda}(\gamma_i, v_{j_k+1})) \quad (31) $$

with $V_{i,j_k,j_l} = (1 - w_{i,j_k})(1 - w_{i,j_l})\delta_{j_k,j_l} + (1 - w_{i,j_k})w_{i,j_l}\delta_{j_k,j_l+1}$ and $W_{i,j_k,j_l} = w_{i,j_k}(1 - w_{i,j_l})\delta_{j_k,j_l+1}$ +
<figure><img src="image_7.png" /><figcaption>Fig. 5: Coronal (top) and axial (bottom) slices of the three CT images reconstructed from the projections used in this study: Forbild phantom (simulated, left) with diaphragm amplitude 20 mm, CIRS phantom (real data, center) and clinical data (real data, right). The two extreme positions are superimposed in green and purple.</figcaption></figure>

$w_{i,j_k} w_{i,j_l} \delta_{j_k+1,j_l+1}$. The Kronecker delta $\delta_{i,j}$ is 1 if $i = j$ and 0 otherwise. These equations assume knowledge of the noise in every pixel of the projections which we estimated with the first order approximation $\text{Var}(g_{\lambda}(\gamma, v)) \approx 1/N(\gamma, v)$ where $N(\gamma, v)$ is the number of photons received in the detector pixel $(\gamma, v)$ [26].

## G. Experiments

1) **Simulated data:** The DCC were validated using simulated acquisitions of the Forbild thorax phantom [27], displayed in Fig. 5, left column. Helical acquisitions were simulated with 360 projections per rotation and four rotations. The duration of a rotation was 0.5 s. The detector size was 920 × 32 pixels of 1.03 × 1.09 mm. The helical pitch was $d = 15.36$ mm (as used in diagnostic 3D CT), the radius of the helical trajectory was 610 mm and the source-to-detector distance was 1113 mm. Ten rays were simulated and averaged per pixel. The acquisitions were centered around $z = 0$ mm. Both noiseless and noisy projections were generated using the reconstruction toolkit (RTK) [28].

To assess the robustness of the proposed metric to noise, several levels of Poisson noise were simulated, from $2.5 \times 10^4$ to $10^5$ photons per pixel in air. This range roughly corresponds to the observed noise level of 3D and 4D patient data. Poisson noise was added to the projections before taking their logarithm. All moments $M_{\tau}$ (15) were computed with $\nu = 0.2$, without adjusting it to the noise level.

When breathing, the phantom's lungs followed the motion $c(t) = A(1 - \cos(2\pi t/T_{\text{resp}}))^2$, with a respiratory period of $T_{\text{resp}} = 2$ s. The breathing started at the beginning of the acquisition and lasted one cycle. The amplitude of the respiratory signal $c(t)$ was also varied in $A \in [0.25, 1]$ mm to assess the sensitivity of the metric to the motion amplitude. In this range, the inferior/superior (IS) peak-to-peak motion was from 5 mm to 20 mm, which covers the average diaphragm motion observed in [29].

2) **Phantom data:** The CIRS Dynamic Thorax Motion Phantom model 008A with a 3 cm diameter tumor in the left lung was scanned with the Somatom go.Sim Siemens

CT scan with a rotation time of 0.35 s and a pitch of 1.34 mm, which is the clinical protocol for 4D CT acquisitions recommended by the manufacturer for the applied tumor motion: a sinusoidal wave of $T_{\text{resp}} = 4$ s with motion amplitudes of 20 mm and 5 mm in the IS and anterior/posterior (AP) directions, respectively. The phantom was scanned twice, once with the tumor static and once with the tumor moving. The 3D CT frames at the two extreme phases of the 4D CT image of the moving phantom are shown in Fig. 5, central column. A reference breathing phase of the phantom data was visually measured by identifying end-exhale projections from the visible displacement of the insert in the projections and computing linear segment from 0% to 100% between two consecutive end-exhale positions. The moments $M_{\tau}$ (15) were computed with $\nu = 0.1$.

3) **Patient data:** One patient, whose end-inhale and end-exhale CT images are shown in Fig. 5, right column, was selected from our database of 4D CT images acquired with the Somatom go.Sim Siemens CT scan. The use of these images respects the French general data protection regulation (GDPR) for medical images. The rotation time was the same as that of the CIRS phantom but the helical pitch was 1.54 mm. The patient projections were truncated because the FOV is not large enough to contain the patient table. Therefore, a pre-processing step consisting in the removal of the table from the projections was necessary to compute the DCC which assume no lateral truncation. We followed the algorithm proposed in [30] which is three steps: first, a 3D reconstruction of the whole volume (body + table) is computed; second, the patient is masked out from the reconstructed image to keep the table part only; third, the table part is forward projected according to the acquisition geometry and the resulting table projections subtracted from the measured projections. To verify that DCC detect consistent projections (assuming a periodic respiratory signal), the reference phase of the patient was visually estimated from the visible displacement of the diaphragm in the projections and was computed in the same way as the phase of the phantom data. The moments $M_{\tau}$ (15) were also computed with $\nu = 0.1$.

## III. RESULTS

### A. Simulated data

1) **Analysis of the DCC:** Fig. 6 illustrates DCC for three pairs of noiseless cone-beam projections of the static acquisition. Corresponding moments perfectly match when the baseline does not cross the FOV (Fig. 6, left). Otherwise, two kinds of situations occur depending on the location of the singular point with coordinates $(\gamma_{\lambda,\lambda}^*, v_{\tau}(\gamma_{\lambda,\lambda}^*))$, i.e. the point of the projection at the intersection with the baseline. This singular point is the same for all B epipolar planes in each projection of the pair (the intersection of the coloured lines at the top of Fig. 6). When the singular points are located in spatially smooth regions in both projections (Fig. 6, middle), the discretization error is small and the corresponding moments match, which
<figure><img src="image_8.png" /><figcaption>Fig. 6: Illustration of the DCC for three pairs of noiseless cone-beam projections (λ, λ'). Top: three pairs of projections and their intersection with five epipolar planes: the plane $Π_{λ,λ'}$ (green), the first and last untruncated planes (blue and purple), and two intermediate planes (red and orange). Bottom: plot of the moments. The columns illustrate DCC when the baseline does not intersect the FOV (left) and when it does (bowtie shape of the colored lines), in a constant region (middle), and close to an edge (right). The colored points match the color of the epipolar plane intersections above.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7: Consistency metric $E_{λ,λ'}$ (top row) and its parts (bottom row) for all pairs $(λ, λ')$ with noisy projections. The blue line and dots correspond to the static acquisition and the orange line and dots correspond to the acquisition where the phantom moves. The dotted black line corresponds to $E_{λ,λ'} = 0.8$. The gray rectangles mark the geometrical situation where the baseline intersects the FOV. This results corresponds to a motion amplitude of 20 mm and the number of photons (per pixel) in air is $10^5$.</figcaption></figure>

indicates that the singularity is correctly dealt with. When it is located around steep spatial gradients of one of the two projections (Fig. 6, right), the error induced by the sampling is significant around the singular point.

Fig. 7 shows the consistency metric $E_{λ,λ'}$ and its different parts for a fixed projection of reference $λ = 4π$ at end-inhale and all other source positions $λ'$ with which DCC can be computed (see Section II-E), for both static and dynamic simulations of the phantom. When the baseline and the FOV do not intersect (white background), the consistency metric for the static acquisition is rather constant and around 0.8. This is the expected value for uncorrelated samples of a normal distribution for which the mean absolute deviation (the numerator of $E_{λ,λ'}$) is $\sigma\sqrt{2\pi^{-1}} \approx 0.8\sigma$ and the standard deviation (the denominator) $\sigma$. The values are centered around a lower value when the baseline and the FOV intersect

(gray background) because the denominator variance (28), calculated under the assumption of uncorrelated moments $M_τ$, is higher than the true variance of the numerator without absolute values due to the correlations of the planes which all intersect in the singular point. When the phantom moves (orange curves with square symbols in Fig. 7), the detection of motion is possible when the baseline does not intersect the FOV (white background) and the inconsistency increases with $|λ - λ'|$ as does the phantom displacement. When the baseline and the FOV intersect (gray background), it is generally difficult to detect motion because the inconsistencies due to motion are dominated by those due to noise and interpolation errors. Motion seems to be detectable only on the side of the gray rectangles, i.e. when the baseline intersects the FOV but not the phantom, in which case the singular point is in an air region of the projection, which is the least noisy and locally constant. For this simulated acquisition with a conventional pitch ($d = 15.36$ mm), the right-hand side of (21) is 7.81 rad. This limit is reached when $|λ' - λ|$ is 4.92 rad, which is slightly above $3π/2 \approx 4.71$ rad (Fig. 7).

2) *Noise and motion sensitivity:* Sixteen acquisitions described in Section II-G.1 were generated for different motion amplitudes and levels of noise (Table I). For each acquisition, the projection of reference $λ = 4π$ was selected, and consistency metrics were computed for all applicable pairs $(λ, λ')$ to obtain a consistency curve as the orange one at the top plot of Fig. 7. To have an overview of the inconsistencies, the consistency metrics are summed up over all applicable pairs of each acquisition in two tables: Table I for the proposed metric $E_{λ,λ'}$ and Table II for a state-of-the-art metric, the mean squared error of the moments over the *B* planes. The means and standard deviations are computed from 10 noise realizations.

For a fixed level of noise, the values of both metrics increased with the motion amplitude. For a fixed motion amplitude, the values of the mean squared error decreased with the noise level (increasing number of photons), as expected. The opposite behavior was observed with the
<table><caption>TABLE I: Sum over all applicable pairs $(\lambda, \lambda')$ for the reference position $\lambda = 4\pi$ of the consistency metric $E_{\lambda,\lambda'}$ for several levels of noise and motion amplitude of the diaphragm A.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Number of photons in air</th></tr><tr><th>2.5 × 10<sup>4</sup></th><th>5 × 10<sup>4</sup></th><th>7.5 × 10<sup>4</sup></th><th>10<sup>5</sup></th></tr></thead><tbody><tr><th>A (mm)</th><td></td><td></td><td></td><td></td></tr><tr><th>0</th><td>269 ± 10</td><td>269 ± 7</td><td>273 ± 5</td><td>282 ± 5</td></tr><tr><th>5</th><td>689 ± 25</td><td>879 ± 21</td><td>1028 ± 15</td><td>1168 ± 24</td></tr><tr><th>10</th><td>1230 ± 28</td><td>1642 ± 24</td><td>1977 ± 23</td><td>2271 ± 27</td></tr><tr><th>20</th><td>2418 ± 15</td><td>3387 ± 25</td><td>4095 ± 33</td><td>4725 ± 25</td></tr></tbody></table>

<table><caption>TABLE II: Sum over all applicable pairs $(\lambda, \lambda')$ for the reference position $\lambda = 4\pi$ of the mean squared error of the moments for several levels of noise and motion amplitude of the diaphragm A.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Number of photons in air</th></tr><tr><th>2.5 × 10<sup>4</sup></th><th>5 × 10<sup>4</sup></th><th>7.5 × 10<sup>4</sup></th><th>×10<sup>5</sup></th></tr></thead><tbody><tr><th>A (mm)</th><td></td><td></td><td></td><td></td></tr><tr><th>0 mm</th><td>353 ± 51</td><td>213 ± 36</td><td>141 ± 26</td><td>120 ± 18</td></tr><tr><th>5 mm</th><td>544 ± 45</td><td>342 ± 37</td><td>303 ± 30</td><td>274 ± 21</td></tr><tr><th>10 mm</th><td>961 ± 81</td><td>771 ± 35</td><td>732 ± 17</td><td>696 ± 39</td></tr><tr><th>20 mm</th><td>2845 ± 158</td><td>2706 ± 89</td><td>2631 ± 56</td><td>2674 ± 42</td></tr></tbody></table>

proposed metric $E_{\lambda,\lambda'}$: for all amplitudes (except when the phantom is static), the values increased with the number of photons. Increasing the number of photons reduces the uncertainty, hence the denominator of $E_{\lambda,\lambda'}$. When the phantom is static ($A = 0$ mm), the values of the proposed metrics are quite constant, unlike the mean squared error which increases with noise. The same applies to the standard deviations when the phantom is moving.

## B. Phantom data

The results of the consistency metric for the CIRS phantom are presented for all applicable pairs with one fixed projection of reference $\lambda = 72\pi$ in Fig. 8 and $\lambda = 100\pi+2.95$ rad in Fig. 9. In Fig. 8, the reference projection was acquired when the tumor was at an extreme position of its motion cycle. In Fig. 9, the reference projection was acquired when the tumor was at an intermediate position of its motion cycle.

For a static tumor (top row in Fig. 8 and Fig. 9), the consistency metric follows the pattern observed in Fig. 7 in the previous section, i.e., rather constant for all pairs with noisier values when the baseline and the FOV intersect. For a moving tumor (middle row in Fig. 8 and Fig. 9), the values of the consistency metric are close to the ones obtained for a static phantom for projections separated by about 4 s for an extreme position (Fig. 8) and about 2 s for an intermediate position (Fig. 9), as expected since the phantom has the same position at mid-exhale and mid-inhale. The reference phase, computed from end-exhale positions of the insert visually identified in the projections, had a period of about 4 s.

Noise and interpolation effects are quite important in both the static and dynamic acquisitions. For pairs whose baselines intersect the FOV, the inconsistencies are mainly caused by the projection noise and those due to tumor

motion are not visible. This is even the case for pairs whose baseline is close to the FOV without intersecting it. This noise can be reduced with a moving average with a $2\pi$ width, using all pairs or only those whose baseline does not intersect the FOV (bottom row in Fig. 8 and Fig. 9). The moving average improves the detection of inconsistencies due to motion in both cases since the noise is averaged out. The values obtained from the static acquisition are constant which confirms the consistency of the projections with their respective reference, while those obtained from the moving acquisitions match the values of the static curves only when the phantom has the same phase as the phase of the reference projection, i.e. when the phantom is at the same position.

The benefit of computing DCC when the baseline intersects the FOV is illustrated by only computing the moving average at points for which the center of the $2\pi$ window has a DCC. Excluding pairs whose baseline does not intersect the FOV results in a dashed orange curve whereas the red curve is more continuous with nearly no gap for one period of motion around the reference value $\lambda$. There are also two more rotations, for $|\lambda - \lambda'| > 25\pi$, for which the only available DCC correspond to pairs with baselines intersecting the FOV (Fig. 8 and Fig. 9).

For this low pitch acquisition ($d = 1.34$ mm), the right-hand side of (21) is 89.63 rad. This limit is first reached when $|\lambda' - \lambda|$ is close to $2\pi$ and is last reached when $|\lambda' - \lambda|$ is 85.43 rad $\approx 27\pi$ which is the range in Fig. 8 and Fig. 9.

## C. Patient data

DCC of patient data are similarly displayed in Fig. 10 for a reference projection at angle $\lambda = 318\pi$ at end-exhale. The respiratory period of the patient was about 3.7 s. The beginning of the phase is taken at the end of expiration. The behavior of the DCC are similar to the real phantom data: the moving average computed from the consistency metric indicates consistent and inconsistent projections that are in agreement with the respiratory phase.

# IV. DISCUSSION

This article derives data consistency conditions for the helical geometry in the detector's coordinates. The change of variables enables the computation of DCC for pairs whose baseline intersects the object (unlike the rebinning approach to a virtual plane [12]), thus increasing the number of pairs used in the consistency metric. This is of significant importance for the detection of local inconsistencies, e.g. isolated blank projections or rapid motion. However, (14) is a Hilbert transform with a singularity when the detector (or equivalently the FOV) and the baseline intersect. This singularity was dealt with by bandlimiting the Hilbert kernel. The results indicate that our implementation detects inconsistencies due to motion when the baseline and the FOV do not intersect (Fig. 6 left and Fig. 7). However, the DCC are sensitive to discretization errors when they do, even without noise (Fig. 6, right), because the Hilbert kernel $h_\nu$ strongly amplifies the
<figure><img src="image_10.png" /><figcaption>Fig. 8: Consistency metric $E_{\lambda,\lambda'}$ computed from a reference projection $\lambda = 72\pi$ acquired when the tumor was at an extreme position of the breathing cycle. Top: static tumor. Center: moving tumor. Bottom: moving average over one gantry rotation of the above plots and tumor motion. The moving average is only computed when the center of the window has a DCC.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9: Consistency metric $E_{\lambda,\lambda'}$ computed from a reference projection $\lambda = 100\pi + 2.95$ rad acquired when the tumor was at an intermediate position of the breathing cycle. Top: static tumor. Center: moving tumor. Bottom: moving average over one gantry rotation of the above plots and tumor motion. The moving average is only computed when the center of the window has a DCC.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10: Consistency metric (top) and moving average (bottom) for the set of patient projections using a reference projection $\lambda = 318\pi$ at end-exhale.</figcaption></figure>

discretization error, despite the Hann windowing, resulting in inconsistencies (Fig. 6, right). This was validated by simulating projections with finer pixels which resulted in lower inconsistencies (e.g. the mean absolute deviation decreases from 1.76 to 0.26 with 10 times smaller pixels for the right pair in Fig. 6). When the baseline and the FOV intersect, DCC are also extremely sensitive to noise. Instead of denoising the projections before computing the DCC, e.g. with an average filter or a median filter, we have multiplied the Hilbert filter with a Hann window to limit the effect of high frequencies, dominated by noise, on the moments. The cut-off frequency $\nu$ should then be adjusted to the noise level of the projections.

Other pairwise DCC [10], [11] based on Grangeat's relation [22] also compare two moments, but these moments do not have a singularity. It was proven that if the projections verify the fan-beam DCC (11), they will also verify the Grangeat DCC, however, the reciprocal is false [31]. The robustness of the different pairwise DCC to discretization and noise might be different since the DCC used in this work have a singularity while the Grangeat DCC require the computation of a derivative (which is also known to be sensitive to noise). Comparing these two DCC was beyond the scope of this work.

Our simulated and real data indicate that DCC can identify breathing motion on a clinical CT scanner. Motion is detected when the baseline does not intersect the FOV but sensitivity to noise might prevent it otherwise (Fig. 7). Sensitivity to noise can be reduced in other pairs with a simple low-pass filter as a moving average (Fig. 8, Fig. 9 and Fig. 10). The proposed metric normalizes the inconsistency by an estimate of the noise level to give more weight to pairs where the motion inconsistency dominates the noise inconsistency in the moving average.

The ability of the consistency metric to detect motion

in real noisy conditions was assessed and we showed that it allowed the detection of small motion, even with a high level of noise (see Table I). The measured inconsistency increased with the amplitude of the motion. Our results on simulated and real data also indicate that motion can be detected in both conventional (3D) and low (4D) pitch acquisitions. However, it is not possible to formally define a minimal amplitude or a velocity at which motion may be detected (for a given pitch and x-ray dose) as the DCC depends on the scanned object or patient. Other sources of inconsistencies, e.g. scatter or beam hardening, might also dominate motion. An evaluation of the derived DCC on an extensive clinical dataset was beyond the scope of this paper.

Due to the large variations of the moments absolute difference with the gantry angle difference of the pair (Fig. 7), it is necessary to normalize this error. The first normalization we considered was by the mean of the two moments. However, pairs with an absolute gantry angle difference close to 180° may have moments around zero since the sign of the integrand may change in the integration interval. This normalization is then inadequate. The impact of noise on the DCC depending on the geometry configuration of the pair led us to normalize the error by an estimate of the noise and we have derived a metric which is adapted to the noise level (constant errors for a static object and constant standard deviation for all the tests in Table I) unlike the mean squared error (increasing errors for a static object and increasing standard deviation when the number of photons decreases for all tests in Table II).

A contribution of this work is the mathematical description of the set of axially-truncated projection pairs of a helical acquisition for which the proposed DCC can be used (II-E). Even if the number of DCC has been increased by dealing with the singularity, axial truncation
of the projections prevents the computation of DCC for some pairs which limits the clinical applications of the DCC in helical CT. A graph approach could be used to calculate the consistency between two projections for which this DCC is not applicable [32]. If one wants to use the DCC to extract a breathing signal, one would need to select a reference projection, e.g. end-exhale to extract all end-exhale positions. But the DCC can be used as is instead, e.g. to sort the projections for 4D CT imaging by maximizing all possible DCC in each subset instead of using an external breathing signal which might not be representative of the internal motion. The DCC could also be used to estimate the parameters of a motion model before reconstruction to compensate for the motion during reconstruction.

V. CONCLUSION

This article derives DCC between pairs of projections acquired on a helical trajectory with a cylindrical detector. The DCC are expressed directly in detector coordinates to compute DCC even if the baseline defined by the two source positions intersects the FOV. In this case, the resulting DCC presents a singularity which we addressed in our implementation. We also identified the set of DCC which can be used with axial truncation. Our tests on simulated and real data demonstrate the ability of the DCC to detect breathing motion.

VI. APPENDIX: MOMENT IN THE CYLINDRICAL DETECTOR COORDINATES

We derive here the expression of the moment $M_τ$ (11) in the detector coordinates $(\gamma, v)$. Given the definition of $c_τ$, which points from the source towards the object, and that we assume no lateral truncation of the projections, the range of $\gamma$ is limited to $[-\gamma_{max}, \gamma_{max}]$. By definition, we have $g_τ(\phi) = g_λ(\gamma, v_τ(\gamma))$. Recalling that the baseline is not orthogonal to the helix axis if the pitch is not zero, i.e. $b_3 \neq 0$ and $c_3 \neq 0$, the denominator of the integrand in (11) is

$$ \begin{aligned} \cos \phi &= \frac{(\mathbf{p}_\lambda(\gamma, v_\tau) - \mathbf{s}_\lambda)}{\sqrt{D^2 + v_\tau^2(\gamma)}} \cdot \mathbf{c}_\tau \quad (32) \\ &= \frac{D}{\sqrt{D^2 + v_\tau^2(\gamma)}} \frac{-b_2 \cos(\gamma - \lambda) - b_1 \sin(\gamma - \lambda)}{n_3} \end{aligned} $$

The last line is obtained by replacing $v_τ$ using (13), the vector coordinates by their values (9) and by using the relation $b_{\lambda, \lambda'} = -c_τ \times n_τ$. The derivative of $\phi$ with respect to $\gamma$ can be calculated by differentiating the expression of $\cos \phi$ with respect to $\gamma$ using the chain rule

$$ \frac{d \cos \phi}{d \gamma} = \frac{d \cos \phi}{d \phi} \frac{d \phi}{d \gamma} = - \sin \phi \frac{d \phi}{d \gamma}. \quad (34) $$

The expression for $\sin \phi$ can be obtained in the same way as $\cos \phi$:

$$ \sin \phi = \frac{(\mathbf{p}_{\lambda}(\gamma, v_{\tau}) - \mathbf{s}_{\lambda})}{\sqrt{D^2 + v_{\tau}^2(\gamma)}} \cdot \mathbf{b}_{\lambda, \lambda'} \quad (35) $$

$$ = \frac{D}{\sqrt{D^2 + v_{\tau}^2(\gamma)}} \frac{c_2 \cos(\gamma - \lambda) + c_1 \sin(\gamma - \lambda)}{n_3}, \quad (36) $$

this time using $c_τ = b_{\lambda, \lambda'} \times n_τ$. Keeping in mind that $v_τ$ is a function of $\gamma$, the differentiation of $\cos \phi$ with respect to $\gamma$ from (33) is tedious but straightforward using the relations $c_τ = b_{\lambda, \lambda'} \times n_τ$ and $b_{\lambda, \lambda'} \cdot n_τ = 0$, resulting in

$$ \frac{d \cos \phi}{d \gamma} = \frac{D^3 (c_2 \cos(\gamma - \lambda) + c_1 \sin(\gamma - \lambda))}{n_3^2 (D^2 + v_\tau^2(\gamma))^{3/2}}. \quad (37) $$

Finally, assembling (34), (36) and (37) yields

$$ \left| \frac{d\phi}{d\gamma} \right| = \frac{D^2}{|n_3|(D^2 + v_{\tau}^2(\gamma))}. \quad (38) $$

The absolute value accounts for the potential change of order of the limits of the integral. Noting $\tilde{g}_λ(\gamma, v_τ(\gamma)) = Dg_λ(\gamma, v_τ(\gamma))/\sqrt{D^2 + v_τ^2(\gamma)}$, the expression of $M_τ$ in detector coordinates is

$$ M_{\tau} = \int_{-\gamma_{\max}}^{\gamma_{\max}} \frac{-\operatorname{sign}(n_3) \tilde{g}_{\lambda}(\gamma, v_{\tau}(\gamma))}{b_2 \cos(\gamma - \lambda) + b_1 \sin(\gamma - \lambda)} d\gamma \quad (39) $$

$$ = \int_{-\gamma_{\max}}^{\gamma_{\max}} \frac{-\tilde{g}_{\lambda}(\gamma, v_{\tau}(\gamma))}{|\cos \alpha| |\cos(\gamma + \bar{\lambda} - \lambda)|} d\gamma \quad (40) $$

where we have injected into the second line the expressions of $n_3$, $b_1$ and $b_2$ in (9) and simplified the resulting expression.

Similarly to (11), (40) might have a singularity if the denominator is null, i.e. if $\gamma = \gamma_{\lambda, \lambda'}^*$ with

$$ \gamma_{\lambda, \lambda'}^* = (\lambda - \bar{\lambda}) \pmod{\pi} - \frac{\pi}{2} \quad (41) $$

$$ = \frac{\lambda - \lambda'}{2} \pmod{\pi} - \frac{\pi}{2}. \quad (42) $$

The derivations of these two expressions of $\gamma_{\lambda, \lambda'}^*$ account for the domain of $|\gamma| < \gamma_{max} < \pi/2$. Geometrically, $\gamma_{\lambda, \lambda'}^*$ is the detector angular coordinate at which the baseline intersects the detector. Injecting (41) into the denominator of (40), we obtain

$$ \begin{aligned} \cos(\gamma + \bar{\lambda} - \lambda) &= \begin{cases} \sin(\gamma - \gamma_{\lambda, \lambda'}^*) & \text{if } (\lambda - \bar{\lambda}) \pmod{2\pi} < \pi \\ -\sin(\gamma - \gamma_{\lambda, \lambda'}^*) & \text{otherwise} \end{cases} \\ &= -\operatorname{sign}(\gamma_{\lambda, \lambda'}^*) \sin(\gamma_{\lambda, \lambda'}^* - \gamma). \end{aligned} \quad (43) $$

The constraint on $(\lambda - \bar{\lambda}) \pmod{2\pi}$ has been translated into sign $\gamma_{\lambda, \lambda'}^*$ by noting that $\gamma_{\lambda, \lambda'}^* > 0$ can be translated into a constraint on $(\lambda' - \lambda)/2$ with (42) and combined with the definition of $\bar{\lambda}$ in (5).
Defining $h(x) = (\pi x)^{-1}$ as the kernel of the Hilbert
transform, we can write

$$
\frac{1}{\sin(\gamma_{\lambda, \lambda'}^* - \gamma)} = \pi h(\sin(\gamma_{\lambda, \lambda'}^* - \gamma)) \quad (44)
$$

$$
= \frac{\pi h(\gamma_{\lambda, \lambda'}^* - \gamma)}{\mathrm{sinc}(\gamma_{\lambda, \lambda'}^* - \gamma)} \quad (45)
$$

where the right hand side is obtained because *h* is homo-
geneous of degree −1. Combining (40), (43) and (45) leads
to the final formulation of the moment given in (14).

ACKNOWLEDGMENT

A CC-BY public copyright license has been applied by
the authors to the present document and will be applied
to all subsequent versions up to the Author Accepted
Manuscript arising from this submission, in accordance
with the grant's open access conditions.

This work was performed using HPC resources from
GENCI-IDRIS.

REFERENCES

[1] H. Barrett and K. Myers, *Foundations of Image Science (1st edition)*. Wiley, 2013.

[2] S. Helgason, *The Radon transform on Euclidean spaces, compact two-point homogeneous spaces and Grassmann manifolds*. The Institute for Advanced Study, Princeton, N. J., USA, 1965.

[3] D. Ludwig, “The Radon transform on Euclidean space,” *Commun. Pure Appl. Math.*, vol. XIX, pp. 49–81, 1966.

[4] Y. Huang, X. Huang, O. Taubmann, Y. Xia, V. Haase, J. Hornegger, G. Lauritsch, and A. Maier, “Restoration of missing data in limited angle tomography based on helgason–ludwig consistency conditions,” *Biomed. Phys. Eng. Express*, vol. 3, no. 3, p. 035015, 2017.

[5] G.-H. Chen and S. Leng, “A new data consistency condition for fan-beam projection data,” *Med. Phys.*, vol. 32, no. 4, pp. 961–967, 2005.

[6] X. Mou, S. Tang, T. Luo, Y. Zhang, and H. Yu, “Beam hardening correction based on HL consistency in polychromatic transmission tomography,” in *Med. Imag. 2008: Phys. Med. Imag.*, vol. 6913, pp. 1026–1035, Int. Soc. Opt. Photon., SPIE, 2008.

[7] T. Würfl, N. Maass, F. Dennerlein, X. Huang, and A. Maier, “Epipolar consistency guided beam hardening reduction-ECC²,” in *Proc. 14th Int. Meeting on Fully Three Dimensional Image Reconstruct. Rad. Nucl. Med.*, pp. 181–185, 2017.

[8] S. Abdurahman, R. Frysch, R. Bismark, S. Melnik, O. Beuning, and G. Rose, “Beam hardening correction using cone beam consistency conditions,” *IEEE Trans. Med. Imag.*, vol. 37, no. 10, pp. 2266–2277, 2018.

[9] S. Tang, K. Huang, Y. Cheng, X. Mou, and X. Tang, “Optimization based beam-hardening correction in ct under data integral invariant constraint,” *Phys. Med. Biol.*, vol. 63, p. 135015, jul 2018.

[10] C. Debbeiler, N. Maass, M. Elter, F. Dennerlein, and T. M. Buzug, “A new CT rawdata redundancy measure applied to automated misalignment correction,” in *Proc. 12th Int. Meeting on Fully Three Dimensional Image Reconstruct. Rad. Nucl. Med.*, pp. 264–267, 2013.

[11] A. Aichert, M. Berger, J. Wang, N. Maass, A. Doerfler, and A. K. Maier, “Epipolar consistency in transmission imaging,” *IEEE Trans. Med. Imag.*, vol. 34, no. 11, pp. 2205–2219, 2015.

[12] J. Lesaint, S. Rit, R. Clackdoyle, and L. Desbat, “Calibration for circular cone-beam CT based on consistency conditions,” *IEEE Trans. Radiat. Plasma Med. Sci.*, vol. 1, no. 6, pp. 517–526, 2017.

[13] Y. Wei, H. Yu, and G. Wang, “Integral invariants for computed tomography,” *IEEE Signal Process. Lett.*, vol. 13, no. 9, pp. 549–552, 2006.

[14] H. Yu and G. Wang, “Motion artifact reduction in fan-beam CT,” *IEEE Trans. Med. Imag.*, vol. 26, no. 2, pp. 249–260, 2007.

[15] M. S. Levine, E. Y. Sidky, and X. Pan, “Consistency conditions for cone-beam CT data acquired with a straight-line source trajectory,” *Tsinghua Sci. Technol.*, vol. 15, no. 1, pp. 56–61, 2010.

[16] R. Clackdoyle, S. Rit, J. Hoskovec, and L. Desbat, “Fanbeam data consistency conditions for applications to motion detection,” in *Proc. 3rd Int. Conf. Image Formation X-Ray Comput. Tomograph.*, pp. 324–328, 2014.

[17] R. Frysch and G. Rose, “Rigid motion compensation in interventional C-arm CT using consistency measure on projection data,” in *Proc. 18th Int. Conf. Med Image Comput. Comput. Assisted interventions*, pp. 298–306, 2015.

[18] R. Clackdoyle, “Necessary and sufficient consistency conditions for fanbeam projections along a line,” *IEEE Trans. Nucl. Sci.*, vol. 60, no. 3, pp. 1560–1569, 2013.

[19] R. Clackdoyle and L. Desbat, “Full data consistency conditions for cone-beam projections with sources on a plane,” *Phys. Med. Biol.*, vol. 58, pp. 8437–8456, 2013.

[20] R. Clackdoyle, L. Desbat, J. Lesaint, and S. Rit, “Data consistency conditions for cone-beam projections on a circular trajectory,” *IEEE Signal Process. Lett.*, vol. 23, no. 12, pp. 1746–1750, 2016.

[21] M. H. Nguyen, *Data Consistency Conditions in 3D Tomography and Scanner Calibration using Analytic Approaches*. PhD thesis, Université Grenoble Alpes, 2021.

[22] P. Grangeat, “Mathematical framework of cone beam 3D reconstruction via the first derivative of the radon transform,” in *Math. Methods Tomograph*. (G. T. Herman, A. K. Louis, and F. Natterer, eds.), pp. 66–97, Springer Berlin Heidelberg, 1991.

[23] S. Patch, “Consistency conditions upon 3D CT data and the wave equation,” *Phys. Med. Biol.*, vol. 47, no. 15, pp. 2637–2650, 2002.

[24] S. Tang, Q. Xu, X. Mou, and X. Tang, “The mathematical equivalence of consistency conditions in the divergent-beam computed tomography,” *Journal of X-Ray Science and Technology*, vol. 20, no. 1, pp. 45–68, 2012.

[25] J. Lesaint, *Data consistency conditions in X-ray transmission imaging and their application to the self-calibration problem*. PhD thesis, Université Grenoble Alpes, 2018.

[26] A. Wunderlich and F. Noo, “Image covariance and lesion detectability in direct fan-beam x-ray computed tomography,” *Phys. Med. Biol.*, vol. 53, no. 10, pp. 2471–2493, 2008.

[27] F. Bergner, T. Berkus, M. Oelhafen, P. Kunz, T. Pan, and M. Kachelriess, “Voxel-based reconstruction combined with motion detection for slow rotating 4D CBCT,” pp. 5107–5112, Oct 2008.

[28] S. Rit, M. V. Oliva, S. Brousmiche, R. Labarbe, D. Sarrut, and G. C. Sharp, “The reconstruction toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the insight toolkit (ITK),” *J. Phys.*, vol. 489, p. 012079, 2014.

[29] S. Rit, M. van Herk, L. Zijp, and J.-J. Sonke, “Quantification of the variability of diaphragm motion and implications for treatment margin construction,” *Int J Radiat Oncol Biol Phys*, vol. 82, pp. e399–e407, Mar 2012.

[30] C. Mory, G. Janssens, and S. Rit, “Motion-aware temporal regularization for improved 4d cone-beam computed tomography,” *Phys. Med. Biol.*, vol. 61, pp. 6856–6877, sep 2016.

[31] J. Lesaint, S. Rit, R. Clackdoyle, and L. Desbat, “GCC and FBCC for linear tomosynthesis,” *Proc. 5th Int. Conf. Image Formation X-Ray Comput. Tomograph.*, pp. 114–118, 2018.

[32] M. Mouchet, S. Rit, and J. Létang, “Motion detection in helical CT using data consistency conditions,” in *2020 IEEE Nucl. Sci. Symp. Med. Imag. Conf. (NSS/MIC)*, pp. 1–3, 2020.