Cone-Beam Pair-Wise Data Consistency Conditions in Helical CT

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.


I. INTRODUCTION
D ATA 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 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], [7], [8], [9], the calibration of geometric parameters such as the detector position and orientation [10], [11], [12], or the detection and the compensation of rigid patient motion [13], [14], [15], [16], [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 Hegason-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 pairwise 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.

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 ortho-normal basis (e x , e y , e z ).The source moves along the helix 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 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 γ ∈ The field-of-view (FOV) is the region defined by The scanner probes the linear attenuation coefficient f of the object in the FOV according to the line integral model 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] and [23].Any pair of source positions (s λ , s λ ′ ) defines such a baseline B λ,λ ′ oriented with the unit vector illustrated in Fig. 2. Fan-beam DCC may be computed in any epipolar plane, i.e., any plane which contains the baseline B λ,λ ′ .The reference epipolar plane λ,λ ′ is defined by its orthogonal Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Fig. 2. 3D illustration of the cone-beam pair-wise geometry.The baseline B λ,λ ′ is defined by the source positions s λ and s λ ′ (with λ < λ ′ in this case).The intersection of the reference epipolar plane Π λ,λ ′ with the two detectors defines the two curves C λ,λ ′ and C λ ′ ,λ (thick dotted lines).The intersection of an epipolar plane Π τ with the two detectors defines the two curves C τ and C τ ′ with τ = (λ, λ ′ , β) and τ ′ = (λ ′ , λ, β) (thick continuous lines).The dotted limits of the detectors are below the reference epipolar plane Π λ,λ ′ and the dashed ones are below the epipolar plane Π τ .
where c λ,λ ′ = cos λ e x + sin λ e y (5) is the unit vector orthogonal to the baseline B λ,λ ′ and the rotation axis e z .The vector c λ,λ ′ points from the midpoint between s λ and s λ ′ (which is on the baseline B λ,λ ′ ) towards the helix axis, with The other epipolar planes τ with τ = (λ, λ ′ , β) are defined by their angle β with the reference plane λ,λ ′ (Fig. 2).Their orthogonal is The angle β is measured positively in the counter-clockwise direction when viewed in the direction b λ,λ ′ .With the third vector ) is an ortho-normal basis which is used in the following to derive cone-beam DCC from fan-beam DCC.If we let ) be the coordinates of these vectors, they can be derived by three successive rotations Fig. 3. Geometric illustration in the epipolar plane Π τ of the angle φ between c τ and the line defined by s λ and p λ (γ, v τ ).Π τ is not orthogonal to the helix axis e z .
with α = arccos(b λ,λ ′ • (e z × c λ,λ ′ )) the angle between the axial planes (defined by z = const) and B λ,λ ′ .We obtain C. Fan-Beam DCC in an Epipolar Plane The intersection of the epipolar plane τ and the cone-beam at source position s λ defines a fan-beam which may be parameterized by the angle φ between c τ and each ray (Fig. 3).We note the line integral at angle φ 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, is constant for all source positions along the line.This integral, which presents two singularities at −π/2 and π/2, must be understood in the sense of the Cauchy principal value.The computation of the moment M τ 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 τ orthogonal to e z by setting Consequently, the two source positions s λ and s λ ′ are the only ones along B λ,λ ′ and ( 11) only provides the following DCC in the epipolar plane τ where M τ ′ is the moment of the projection g τ ′ with τ ′ = (λ ′ , λ, β) and β in a range derived in (II-E).
We derive the expression of the moment M τ in the detector coordinates (γ , v).At source position s λ , the epipolar plane intersects the detector along the curve C τ defined by where ).The ray direction p λ (γ , v τ (γ )) − s λ now depends on γ only.We have a one-to-one mapping between the angles φ and γ and we can change variable from φ to γ 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 where is the Hilbert kernel and γ * λ,λ ′ is a singularity corresponding to the angular coordinate at which the baseline may intersect the detector.

D. Discretization
The detector has N cols × N rows pixels of size γ × v .Their coordinates γ and v are 14) is discretized as follows: where The projection value is obtained with linear interpolation with 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 where ν ∈ (0, 1] is a fraction of Nyquist's frequency and controls the bandwidth of the kernel.

E. Applicability Intervals for α and β
The DCC are only valid if the fan-beam projections g τ (φ) (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 λ − λ ′ ) and some epipolar planes (i.e.some β) for which the detector curves C τ and C τ ′ are between the inferior and superior limits of the detector (see Fig. 2): If this inequality is not verified for one of the two curves C τ or C τ ′ , as illustrated in Fig. 4, the resulting fan-beam projections g τ and g τ ′ are truncated and the DCC cannot be used.Replacing (n 1 , n 2 , n 3 ) in v τ (γ ) (13) by their expression in (9) and λ − λ by γ * λ,λ ′ from (41), we have, similarly to (43), sin(γ + λ − λ) = − sign γ * λ,λ ′ cos(γ − γ * λ,λ ′ ) and the previous inequality becomes (19) which gives the following conditions on the angles α and β: This constraint on α indirectly constrains the interval separating two source positions where DCC can be applied because α increases with the angular distance |λ ′ − λ| between the pair of source positions.This specific condition where γ coincides with the singularity, i.e. the direction of the baseline, does not depend on β since β is the angle defining the epipolar planes around the baseline. Using which shows that the smaller the helical pitch d, the larger the length of the applicability interval.The sine function in this Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
inequality changes the applicability interval from rotation to rotation.Almost full rotations are applicable ((21) verified) for the rotations around the reference position λ (i.e. when |λ ′ −λ| 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 (λ ′ − λ) mod 2π = π with a length decreasing when the angular distance to the reference where and are the equations of the smallest and highest angle β for which the epipolar curve v τ (γ ) is not truncated for each detector column γ .Because γ * λ,λ ′ is not the same for s λ and s λ ′ , ( 22) must be verified for the two source positions.However, it can be shown that l λ ′ ,λ (γ ) = −u λ,λ ′ (−γ ) and u λ ′ ,λ (γ ) = −l λ,λ ′ (−γ ) because a π rotation around the axis defined by (s λ + s λ ′ )/2 and c λ,λ ′ swaps the two detectors.We therefore obtain the condition: If β max ≤ 0, the detector is too small to compute a DCC for this pair of source positions.Otherwise, we sample β in B = ⌊2β max D/ v ⌋ planes.

F. Consistency Metric
We define the consistency metric between two cone-beam projections as The numerator of E λ,λ ′ 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 M λ,λ ′ and M λ ′ ,λ 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 The √ B weight in the denominator of E λ,λ ′ makes it independent of B, i.e., compensates the decrease of Var(M λ,λ ′ ) 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 τ depends on the weight of each projection value g λ (γ , v) which also strongly varies with τ and γ in the computation of the moments M τ in (15).
The variance Var(M λ,λ ′ ) 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 The terms in the sums above are expanded using (15).We note and obtain with The Kronecker 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 Var (g λ (γ , v)) ≈ 1/N (γ , v) where N (γ , v) is the number of photons received in the detector pixel (γ , 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.09mm.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 × 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 τ (15) were computed with ν = 0.2, without adjusting it to the noise level.
When breathing, the phantom's lungs followed the motion c(t) = A(1 − cos(2πt/T resp ) 2 ), with a respiratory period of T 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 ∈ [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 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 τ (15) were computed with ν = 0.1.
3) Patient Data: One patient, whose end-inhale and endexhale 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 τ (15) were also computed with ν = 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 (γ * λ,λ ′ , v τ (γ * λ,λ ′ )), 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 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 σ √ 2π −1 ≃ 0.8σ and the standard deviation (the denominator) σ .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 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.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 .
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 ≈ 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 proposed metric E λ,λ ′ : 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 λ,λ ′ .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 λ = 72π in Fig. 8 and λ = 100π + 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 Fig. 10.Consistency metric (top) and moving average (bottom) for the set of patient projections using a reference projection λ = 318π at end-exhale.
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π 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π 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 λ.There are also two more rotations, for |λ − λ ′ | > 25π, 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 |λ ′ − λ| is close to 2π and is last reached when |λ ′ − λ| is 85.43 rad ≈ 27π 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 λ = 318π 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 ν strongly amplifies the 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 ν 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 depend 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 (γ , 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 γ is limited to [−γ max , γ max ].By definition, we have g τ (φ) = g λ (γ , v τ (γ )).Recalling that the baseline is not orthogonal to the helix axis if the pitch is not zero, i.e. b 3 ̸ = 0 and c 3 ̸ = 0, the denominator of the integrand in (11) is cos φ = ( p λ (γ , v τ ) − s λ ) ( The last line is obtained by replacing v τ using (13), the vector coordinates by their values The expression for sin φ can be obtained in the same way as cos φ: The constraint on (λ − λ) (mod 2π) has been translated into sign γ * λ,λ ′ by noting that γ * λ,λ ′ > 0 can be translated into a constraint on (λ ′ − λ)/2 with (42) and combined with the definition of λ in (5).
Defining h(x) = (π x) −1 as the kernel of the Hilbert transform, we can write where the right hand side is obtained because h is homogeneous of degree −1.Combining (40), ( 43) and (45) leads to the final formulation of the moment given in ( 14).

Fig. 1 .
Fig. 1.Cone-beam geometry of the helical scanner used in this study.

Fig. 4 .
Fig.4.Example of a plane Π τ for which the DCC cannot be computed because the curves C τ and C τ ′ do not verify the condition(18).The truncated part of the curves are shown in red.

Fig. 5 .
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.

Fig. 6 .
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.

Fig. 8 .
Fig. 8. Consistency metric E λ,λ ′ computed from a reference projection λ = 72π 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.

Fig. 9 .
Fig. 9. Consistency metric E λ,λ ′ computed from a reference projection λ = 100π + 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.

( 9 )
and by using the relation b λ,λ ′ = −c τ × n τ .The derivative of φ with respect to γ can be calculated by differentiating the expression of cos φ with respect to γ using the chain rule d cos φ dγ =

TABLE I SUM
OVER ALL APPLICABLE PAIRS (λ, λ ′ ) FOR THE REFERENCE POSITION λ = 4π OF THE CONSISTENCY METRIC E λ,λ ′ FOR SEVERAL LEVELS OF NOISE AND MOTION AMPLITUDE OF THE DIAPHRAGM A TABLE II SUM OVER ALL APPLICABLE PAIRS (λ, λ ′ ) FOR THE REFERENCE POSITION λ = 4π OF THE MEAN SQUARED ERROR OF THE MOMENTS FOR SEVERAL LEVELS OF NOISE AND MOTION AMPLITUDE OF THE DIAPHRAGM A