Exact computation of projected sphere centres in cone beam x-ray projections

A method is presented for fitting the projected centres of spheres in cone beam x-ray imaging. By using a suitable coordinate system, the method allows direct and exact calculation of the sphere centre without fitting the projection shape with an ellipse and correcting from the ellipse centre to the sphere centre. Advantages in numerical implementation result from the number of unknown variables being reduced compared to ellipse fits. Additionally, the orientation of the detector relative to the x-ray source can be obtained from fitting the shapes of projections of multiple spheres without knowledge of the positions or dimensions of the spheres. The accuracy of the method is compared to other techniques using simulated x-ray projections.


Introduction
X-ray computed tomography (CT) has become a valuable tool for metrology and non-destructive testing of mechanical parts with complex geometries and is used in a wide range of industrial applications.It allows capturing an entire object in one scan including internal structures, which may be inaccessible to tactile measurement devices or optical instruments using visible light.However, reliable quantitative measurements from CT scans require an accurate calibration of the geometry of the CT system [1].
For cone beam x-ray CT systems, a common geometry calibration method is to analyse projection images obtained from a calibrated reference object consisting of a set of markers [2,3]: by measuring the positions and dimensions of the projected marker objects in the x-ray projection images, the projection Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.magnification can be inferred, which ultimately links image pixels to actual object dimensions.
Due to their symmetry, spheres are a natural choice for marker objects with wide application as components of x-ray CT reference objects.The projection of a sphere in a cone beam x-ray CT system with a flat detector has an elliptical shape [4] and can easily be fitted.However, it has the complicating property that the centre of the projection ellipse does not coincide with the projection of the sphere centre [4], which represents the position of the sphere.For accurate geometry calibration using spheres as reference object markers, this deviation between the projection ellipse centre and the projected sphere centre needs to be corrected in general.A first correction was derived by Deng et al [5] to infer the projected sphere centre from the ellipsoid centre.However, this correction is numerically unstable due to a division by zero close to the detector centre and therefore cannot be universally applied without precaution.An adapted correction formula was proposed by Butzhammer et al [6] based on an approximation that eliminates the numerical instability of the Deng method.
In this work, we present a different approach that does not rely on fitting ellipses to the sphere projection image but transforms the detector image into a suitable coordinate system where the projection can be fitted with an equation describing the edge of the sphere projection.The sphere centre projection can then be calculated directly from a set of fit parameters, requiring only the distance between the x-ray source and the detector (sdd) as a supplementary input quantity.A novelty of our method is that in addition to accurately providing the position of the sphere centre projection, it can be used for measuring the orientation of the detector plane with respect to the x-ray source by evaluating multiple sphere projections with different sphere positions but constant source-detector configuration.

Geometry description in Euclidian space
For the description of the sphere projection geometry, we define a coordinate system with an ideal (point-like) x-ray source at the origin and a detector in the plane perpendicular to the x-axis at coordinate x = sdd (source detector distance).An ideal sphere of radius r is placed at coordinates as depicted in figure 1. Due to the symmetry of a sphere, the projection of the sphere edge must be described by a cone that is tangential to the sphere and has its apex at the origin (x-ray source).The intersection of this tangent cone with the sphere edge must be a circle and all points on this circle have the same distance r t = √ R 2 0 − r 2 from the origin.The points on the tangent circle must therefore satisfy 2 equations for the distance from the sphere centre and from the origin: Equation ( 2) can be rewritten as (3)

Definition of 'directional coordinates'
For bringing the tangent circle equation into a usable form, we use a coordinate system that we call directional coordinates and define as where the definition of the angles φ , θ shown in figure 2 is not necessary but illustrates the link to the more commonly used spherical coordinates.The purpose of this choice is that any point of the sphere and its projection onto the detector have the same η, ζ coordinates in this coordinate system.An illustration of the conversion of a circle and an ellipse in Cartesian coordinates to directional coordinates is shown in figure 3.
Cartesian coordinates are expressed in directional coordinates as And therefore

Tangent circle equation in directional coordinates
Expressing equation ( 1) in directional coordinates and applying equation (3) yields Using the fact that R = r t for all points of the tangent circle and dividing by r 2 t , equation ( 7) becomes Equation ( 8) describes in directional coordinates both the cone tangent circle of a sphere and the edge of the sphere's projection on the detector with 3 fit parameters These parameters do not yet correspond to the sphere centre in directional coordinates, but they can be converted using the relation The sphere centre in directional coordinates and the sphere radius relative to the sphere distance from the origin are therefore obtained as From the directional coordinates of the sphere centre, the Cartesian coordinates [x p , y p , z p ] of the sphere centre projection can be calculated from equations ( 5) and ( 6) as ) 2 (11) and analogously, the sphere radius magnified to sdd is For fitting the projection of a sphere with equation (8) in directional coordinates using the coordinate transformation in equation ( 4), the detector pixel grid first has to be converted into Cartesian coordinates as described in the following section.

Conversion pixels ⇔ Cartesian coordinates
Note that while a general fit of an ellipse (as the general shape of a sphere projection) in a detector image uses 5 fit parameters (centre y and z coordinates, major and minor axes, orientation angle), equation (8) appears to have only 3 degrees of freedom.The 2 missing parameters are 'hidden' in the coordinate system and given by the pixel coordinates [u 0 , v 0 ] where the xray beam is perpendicular to the detector plane, defined as the intersection of the x-axis with the detector plane in Cartesian coordinates (y = z = 0) in figure 1.In a conventional, wellaligned x-ray imaging setup, this perpendicular pixel will normally be close to the detector centre, but it can deviate from the centre in general.
Using the pixel size (side length) d px for a given detector and the sign convention defined in figure 2, the conversion from image pixel coordinates [u, v] to Cartesian coordinates of the imaging setup (as defined in figure 1) is given by Conversely, the pixel coordinates [u p , v p ] of the projected sphere centre are obtained as For an x-ray projection image containing n (spatially separated) spheres, fitting all sphere projections with ellipses requires 5n fit parameters (or 4n + 2 if the orientation of the ellipses is restricted to point towards the image centre).The method described here, however, only requires 3n + 2 parameters, since the orientation of the detector with respect to the x-ray beam is independent of the number or position of spheres that are placed in the cone beam and the eccentricity of the projection ellipse is incorporated in the coordinate system.
Using a projection of multiple spheres (or multiple projections of variable sphere positions with fixed detector orientation) can therefore allow for fitting the perpendicular pixel coordinates [u 0 , v 0 ] from the shapes of the projected spheres without knowing the actual sphere positions or diameters.

Sphere edge parametrization
In some cases, it may be more convenient to model the position of the edge of a sphere projection rather than its centre.For instance, if the projection of a sphere is overlapped by some sort of support structure to one side of the sphere, fitting the edge pixels on the support side will be inaccurate or unstable and influence the fitted sphere centre, but the edge pixels on the clear side may still be fitted accurately.
In order to be able to model the edge of a sphere projection with subpixel resolution, a parametrization of the tangent circle of the sphere was derived.
Using geometric similarity (or trigonometric identities and the half cone angle θ r ) in figure 4, it can be seen that And using equations ( 9) and ( 10) The unit vector ⇀ e 0 pointing from the origin to the sphere centre ⇀ R 0 can be expressed in directional coordinates of the sphere centre: An arbitrary unit vector With Since ⇀ e ⊥ was defined with a normalization factor, we can choose e ⊥,y = 1 and obtain This orthogonal vector ⇀ e ⊥ can be rotated around ⇀ e 0 by an angle α using Rodrigues' rotation formula: The parametrization of the tangent sphere edge in Euclidian space is then given by And in directional coordinates All quantities required to calculate this expression (and transform to pixel coordinates) can be obtained from the sphere projection fit and allow easy numerical calculation of for instance the [u, v] extrema of the sphere edge protection with interpolated (fitted) sub-pixel resolution.

Application example
To illustrate the application of the fit method, a projection of a steel sphere is shown in figure 5.The edge of the projected sphere is obtained by calculating the absolute value of the image gradient and binarizing the gradient using a threshold value, resulting in a discrete set of edge pixel coordinates.These pixel coordinates are transformed into directional coordinates and fitted with equation (8).

Accuracy evaluation of the method
In the following, all projection shape fits were performed by fitting the shape edge with geometrical objects using a simple image gradient threshold for edge detection.
For assessing the accuracy of the directional coordinates fit, projections of a regular array of spheres in an 8 mm grid were simulated using aRTist [7] software (figure 6).The central sphere was slightly decentred in order to avoid division by zero but remain in the range of numerical instability for the Deng correction on purpose.The detector was modelled with 4000 × 4000 px with 0.1 mm pixel size.Spheres were simulated with diameters of 2 and 4 mm.The source detector distance (sdd) was varied between 500 and 1500 mm, with the simulation magnification kept constant at 10 by adjusting the source object distance accordingly.
The inaccuracy of the ellipse fit can be on the order of a few pixels for a large cone beam angle (short sdd), whereas the  Figure 7 shows a comparison of directional coordinates and corrected ellipse fits for different numerical simulation parameters.It can be seen that the influence of the simulation settings on the sphere centre accuracy is larger than the fit/correction method, except for the central sphere where the Deng correction suffers from numerical instability.
Interestingly, adding image noise appears to improve the accuracy of the fits.A possible explanation for this observation  might be that the fit accuracy is limited by the discrete pixel resolution and adding noise to it introduces an anti-aliasing effect.
The mean absolute centre deviation for all spheres per projection image is summarized for noise-free simulations in figure 8. None of the methods can be called systematically best, but the Deng correction appears to perform slightly worse than directional coordinates and the Butzhammer correction on average (note that the numerical instability of the Deng correction for the central sphere is reflected in the average deviation).
The dependence of the sphere fit accuracy on the cone beam angle is shown in figure 9.A slight increase of the fit deviation can be seen at large cone beam angles for the smaller sphere diameter, but the influence is negligible over a wide range of angles.

Application for determining the CT detector orientation
The evaluation method was applied to the characterisation of the orientation of the detector in a cone beam CT system by fitting the perpendicular pixel coordinates [u 0 , v 0 ] from both simulated and measured sphere projection images.

Perpendicular pixel fit from simulated data
For testing the accuracy of fitting the perpendicular pixel coordinates with an exactly known geometry, a projection of an irregular array of spheres was simulated as shown in figure 10 for different sdd values.To prevent a bias due to a horizontally asymmetric sphere arrangement, the sphere array was horizontally mirrored for a second projection (corresponding to a 180 • rotation around a vertical axis in a CT setup).
In order to introduce a horizontal drift in the perpendicular pixel coordinates with changing source detector distance, the detector was rotated around a vertical axis by 0.1 • for all simulations.
The result of the perpendicular pixel coordinates fit is shown in figure 11 without image noise and in figure 12 with simulated image noise.The accuracy of the fit is approximately on the order of one pixel, with decreasing precision as the source detector distance increases and the cone beam angle  decreases.Adding image noise has an insignificant effect on the fit.

Perpendicular pixel fit from measurements
For testing the perpendicular pixel fit on real data, a number of Ø 1 mm steel spheres were glued on a 2 mm thick carbon fibre plate as shown in figure 13, where neither the sphere arrange-ment nor the spheres themselves were calibrated.Again, two projections offset by a 180 • rotation were recorded for multiple sdd values.
The result of the perpendicular pixel fit is shown in figure 14.From the slope of the linear fits, a mean detector angle of 0.06 • (horizontal) and −0.02 • (vertical) relative to the detector movement axis can be calculated.

Conclusion
The calculation method for the sphere centre projection coordinates presented in this work has been verified by comparison to state-of-the-art approximate methods.The difference between the method reported here and the Butzhammer correction is smaller than the numerical accuracy of projection simulation and shape fitting.
From projections of multiple spheres, the pixel coordinates where the x-ray beam hits the detector in a right angle can be determined with an accuracy that can reach the order of the detector pixel size, depending on the number of spheres placed in a projection and on the cone beam angle.The method therefore allows obtaining information about the geometry of an xray imaging setup without requirement for calibrated reference objects.

Data availability statement
The data cannot be made publicly available upon publication because no suitable repository exists for hosting data in this

Figure 1 .
Figure 1.Definition of the sphere projection geometry.

Figure 2 .
Figure 2. Definition of the image pixel coordinate system [u, v] and illustration of the angles corresponding to spherical coordinates.

Figure 3 .
Figure 3. Conversion of a circle and an ellipse from Cartesian (left) to directional (right) coordinates.Note that the blue circle is centred at (y,z) = (0,0) and therefore is circular in both coordinate systems.

Figure 4 .
Figure 4. Illustration of the rotational parametrization of the tangent circle of the sphere edge and its projection cone.

Figure 5 .
Figure 5. Projection image of a steel sphere (a), image gradient (b), binarized gradient (c) and edge fit in directional coordinates (d).

Figure 6 .
Figure 6.Example projection of the Ø 4 mm sphere array with source-object distance 50 mm and sdd 500 mm.

Figure 7 .
Figure 7.Comparison of simulated sphere centre fit accuracy for directional coordinate fit and corrected ellipse fit methods at sdd = 500 mm with different projection simulation settings: no image noise (left), exaggerated image noise (centre) and no multisampling (no subpixel simulation resolution, right).

Figure 8 .
Figure 8. Mean absolute sphere centre deviation with the standard deviation as error bars for different correction methods and directional coordinates for simulated spheres of 2 mm (left) and 4 mm (right) diameter.

Figure 9 .
Figure 9. Absolute sphere centre deviation as a function of the cone beam angle for different correction methods and directional coordinates for simulated spheres of 2 mm (left) and 4 mm (right) diameter.

Figure 10 .
Figure 10.Example projection of the simulated irregular sphere array used for fitting the perpendicular detector pixel coordinate.

Figure 11 .
Figure 11.Fitted perpendicular pixel coordinates from simulated projection images without noise in horizontal and vertical offset from the image centre.

Figure 12 .
Figure 12.Fitted perpendicular pixel coordinates from simulated projection images with noise in horizontal and vertical offset from the image centre.

Figure 13 .
Figure 13. 2 mm thick carbon fibre plate with Ø 1 mm steel spheres glued on top in an irregular pattern.

Figure 14 .
Figure 14.Horizontal and vertical deviation of the perpendicular pixel coordinate from the image centre as a function of the source detector distance extracted from measured x-ray projection images.