Offline generator for digitally reconstructed radiographs of a commercial stereoscopic radiotherapy image‐guidance system

Abstract Purpose Image‐guided radiotherapy (IGRT) research sometimes involves simulated changes to patient positioning using retrospectively collected clinical data. For example, researchers may simulate patient misalignments to develop error detection algorithms or positioning optimization algorithms. The Brainlab ExacTrac system can be used to retrospectively “replay” simulated alignment scenarios but does not allow export of digitally reconstructed radiographs (DRRs) with simulated positioning variations for further analysis. Here we describe methods to overcome this limitation and replicate ExacTrac system DRRs by using projective geometry parameters contained in the ExacTrac configuration files saved for every imaged subject. Methods Two ExacTrac DRR generators were implemented, one with custom MATLAB software based on first principles, and the other using libraries from the Insight Segmentation and Registration Toolkit (ITK). A description of perspective projections for DRR rendering applications is included, with emphasis on linear operators in real projective space P3. We provide a general methodology for the extraction of relevant geometric values needed to replicate ExacTrac DRRs. Our generators were tested on phantom and patient images, both acquired in a known treatment position. We demonstrate the validity of our methods by comparing our generated DRRs to reference DRRs produced by the ExacTrac system during a treatment workflow using a manual landmark analysis as well as rigid registration with the elastix software package. Results Manual landmarks selected between the corresponding DRR generators across patient and phantom images have an average displacement of 1.15 mm. For elastix image registrations, we found that absolute value vertical and horizontal translations were 0.18 and 0.35 mm on average, respectively. Rigid rotations were within 0.002 degrees. Conclusion Custom and ITK‐based algorithms successfully reproduce ExacTrac DRRs and have the distinctive advantage of incorporating any desired 6D couch position. An open‐source repository is provided separately for users to implement in IGRT patient positioning research.


INTRODUCTION
The ExacTrac system (Brainlab AG, Munich, Germany) is an image-guidance system designed for spatially precise image-guided radiotherapy (IGRT) treatments such as stereotactic body radiotherapy (SBRT) and intracranial stereotactic radiosurgery. A detailed description of the ExacTrac system may be found elsewhere. 1 A brief overview is as follows. The X-ray system consists of two X-ray imagers with floor-mounted kV sources and ceiling-mounted stereoscopic flat-panel detectors. System geometry is illustrated in Figure 1. Each imager has a field of view (FOV) of approximately 10 cm at isocenter. The system can be used to determine a shift with six degrees of freedom (6D), that is, a threedimensional translation and pitch, yaw, and roll rotation to bring the radiotherapy target position at the time of treatment as close as possible to the planned position in the coordinate system of the treatment machine. The ExacTrac system is typically installed in conjunction with C-arm gantry radiotherapy machines, although other configurations are possible. 2,3 Since the imaging source and detector are mounted to the floor and ceiling, respectively, they are not directly attached to the C-arm. Hence, there is no geometric calibration required for gravity-induced mechanical flex, which is an advantea of ExacTrac over traditional C-arm imagers. 4 Translational and rotational shifts are optimized by iteratively minimizing a gradient-based cross-correlation similarity statistic between the X-ray image and dynamically generated digitally reconstructed radiographs (DRRs) F I G U R E 1 Linac room geometry with ExacTrac image-guided radiotherapy (IGRT) system installed. X-ray tube 1 is paired with flat panel detector 1, and likewise tube 2 is paired with detector 2. The central beamlines intersect at the isocenter at a crossing angle φ. The oblique plane spanned by the central beamlines is at an angle of incline θ with respect to the floor derived from the planning computed tomography (CT) scan. 5 The X-ray system is coupled to an infrared camera system that is used for precise positioning of a custom carbon-fiber treatment couch. At our institution, the ExacTrac system is mounted with Varian TrueBeam, TrueBeam STx, and Novalis Tx radiotherapy devices and is used predominantly in the treatment of targets that are rigidly registered to bony anatomy such as intracranial, spinal, and bone metastasis targets.
In order to study the behavior of DRRs with arbitrary 6D CT transformations, it may be necessary to perform offline DRR rendering. For example, at our institution, we are developing algorithms to detect patient positioning errors in IGRT images. 6,7 Algorithm training required the use of simulated errors produced by matching clinically acquired X-ray projections with DRRs of adjacent ("wrong") vertebral bodies. The ExacTrac offline preparation and review station allows simulation of alignment to the wrong vertebral body but does not allow for export of the corresponding optimized DRRs. Thus, our research required generating DRRs that would accurately mimic those produced by the ExacTrac system for arbitrary patient positions. Here, we present two DRR generators built with custom software, one programmed from scratch and the other incorporating open-source image processing C++ libraries from the Insight Segmentation and Registration Toolkit (ITK). 8 The user may input either assumed in-room measurements or exact configuration values into the generators. We will describe how to obtain ExacTrac renderer values and convert them to suitable parameters for DRR generation with ray tracing. Finally, we will validate our algorithms by comparing our generated DRRs to reference DRRs from the ExacTrac system using landmark analysis and rigid registrations to identify any systematic translational shifts. DRRs of objects in treatment position are exportable from the ExacTrac system and serve as a reference in assessing whether our methods successfully reproduce the ExacTrac geometry. By contrast, the ExacTrac system will not export DRRs with couch shifts that were not actually used at the time of treatment. Our research overcomes this limitation. We believe the methodology presented below to simulate the ExacTrac system geometry, along with corresponding open-source software released by our group, may be useful to other groups engaged in IGRT research. Our software may be found at https://github.com/jchartersmp/DRR.

METHODS
Typically, the CT dataset is saved in digital imaging and communications in medicine (DICOM) format. In our calculations, it is required to convert between machine-centric and patient-centric coordinate systems described as follows. The machine-centric coordinates are given by the International Electrotechnical Commission (IEC), where +X is right, +Y is forward, and +Z is up.The second coordinate system is patient-centric,with our scans specifically in a head-first supine (HFS) order, where +X is right, +Y is down, and +Z is forward. If a subject is scanned in a different position, the user may need an additional coordinate conversion. The Exac-Trac renderer values obtained below assume IEC coordinates,whereas the DICOM headers,and consequently the ITK-imported dataset, assume patient-centric coordinates. We refer to X-ray tube 1 as the tube on the −X side of the linac, which is paired with flat panel detector 1 on the +X side of the linac (in either IEC or HFS coordinates). Tube 2 and detector 2 are on the opposite sides accordingly.

DRR generator from first principles
Given a CT volume dataset and an RT Plan DICOM, it is possible to create stereoscopic DRRs with only four in-room measurements. The necessary values are source-to-image distance (SID), source-to-object distance (SOD), oblique plane angle , and central beamline crossing angle . By symmetry of the ExacTrac system, it often suffices to consider one tube-detector pair, say tube 1 and detector 1. The SID is then the distance from tube 1 to the center of detector 1. The SOD is the distance from tube 1 to the linac isocenter. Each X-ray beam originates at a point source in its X-ray tube and proceeds diagonally upward to its corresponding flat panel detector. We suppose the beams are composed of beamlets, or individual rays, that are in bijection with the individual DRR pixels. We refer to the ray that strikes the center of its corresponding flat panel as the central beamline. The angle between the central beamlines is referred to as the crossing angle . The oblique plane is defined as the unique plane containing the pair of central beamlines. We define the oblique angle as the angle of incline between the oblique plane and the floor.
The central beamline 1 direction d 1 is given by which is then normalized for convenience. Here, R (X ) is a counterclockwise rotation matrix about the IEC Xcoordinate. Since a unit direction effectively parameterizes the beamlet by path length, the tube 1 focal point t 1 is found by solving which is substituted into to get the flat panel 1 center c 1 . By symmetry, the central beamline 2 direction d 2 and tube 2 position t 2 are easily found by sending X ↦ −X . Assume that each central beamline is orthogonal to its corresponding flat panel. Now we need to find the orthogonal in-plane flat panel x and y directions. The procedure is identical for both tube-detector pairs. The first step is to project the central beamline directions onto the floor that is, set the Z-components to zero, and renormalize.
Here, the subscript indicates that the equation holds for either tube-detector pair. The normalized floor vectors d ′ 1,2 are then rotated about the Z axis The flat panel y directions p Finally, the in-plane x directions p to yield right-handed detector coordinate systems. The X-ray tube focal points, detector positions, and detector axes are brought into a custom-made raytrace DRR generator. We implemented the Siddon-Jacobs ray tracing algorithm to calculate radiological paths incident on the detectors. 9,10 All scripts are written in the MATLAB environment (MathWorks, Natick, MA). We assumed a linear correspondence between CT number and voxel intensity. The actual mapping displayed on the ExacTrac workstation is proprietary to Brainlab.
For verification,a simple DRR generator using in-room measurements is presented. Using the oblique and centerline angles in, 11 together with the SID and SOD from the ExacTrac clinical user guide, we reproduced nearly perfect DRRs for our phantom.

DRR generator with ExacTrac configuration and ITK
ExacTrac stores information on all imaged subjects, including the CT volume, DRRs for each fraction, treatment logs indicating X-ray correction shifts, as well as system geometric information. In particular, the information for perspective projection transformations is contained in the configuration (initialization) file. We confirmed our methodology on ExacTrac versions 6.1.1, 6.2.1, and 6.2.3 (in use from 2014 to present at our institution). Within the configuration file, there exists a section "FlatPanel" containing keys "MLinToFlat1" and "MLinToFlat2" for flat panels 1 and 2, respectively. The value of each key is a comma-separated string of numbers. The first number is always a zero, which should be disregarded. In order to convert the remaining numbers into a linear transformation, we rearrange them into a 3 × 4 matrix. Row-major order is required, so that each row is completely filled before moving to the next row. Henceforth, we refer to the constructed matrix as the renderer matrix, as it encapsulates all the information for DRR rendering. The renderer matrices for both detectors are expressed in IEC coordinates.
Given CT slices stored as DICOM files, we use a built-in ITK function that imports the CT dataset into a nearly raw raster data file format. 8 After importing into ITK, the CT volume must be translated by subtracting out the isocenter, which is found in the RT Plan DICOM. This translation effectively moves the CT isocenter to the linac isocenter, which corresponds to the IEC and HFS origin. Suppose that a DRR render window is initialized on the origin, facing the HFS Z direction. We seek linear transformations that convert the camera to the X-ray focal points for 2D flat panel viewing.
An overview of projective geometry is described in Angel. 12 It is helpful to regard linear operators as 4 × 4 matrices on real projective space ℙ 3 . Given a vector in Euclidean space ℝ 3 , we include a homogeneous coordinate to replace it with a vector in ℙ 3 . We will describe all the transformations that constitute a renderer matrix M given by The first matrix is a translation transformation T,which sends the camera from the origin to the X-ray tube focal point. We have where t (i) is the i th IEC component of the focal point position t. In order to aim the camera at the flat panels, we then need to apply a rotation operator R, which is a homogenized version of the flat panel transposed direction cosines. We have where p (j)(i) is the i th IEC component of the j th in-plane flat panel direction.
Once the camera translations and rotations are applied,we may express a projectional calibration matrix P in its standard form A scaling matrix S det is required to convert camera coordinates into detector coordinates. In our ExacTrac system, the detectors have a standard 200 × 200 mm 2 FOV, defined at the flat panels, and 512 × 512 pix 2 resolution after imposing 2 × 2 binning. It follows that the pixel spacings s x and s y are 0.39 mm∕pix. If the data do not undergo 2 × 2 binning, then the pixel spacings must be divided by two accordingly. We also have pixel shifts c x and c y that represent where the upper left corner of the DRR is located with respect to the detector center. The general form of a scaling matrix is Finally, it is often convenient to convert to normalized device coordinates (NDC) in the renderer, where a frustum is mapped to the unit cube. In DRR rendering, we impose a simplified NDC scaling matrix S NDC for a frustum that depends on the detector FOV and SID. The desired matrix is where detector width W and height H are expressed in pixels. When both scaling matrices are applied to the standard calibration matrix, we arrive at Notice that the third and fourth rows are identical because depth information is lost in 2D projections. Thus, the fourth row may be removed to yield a 3 × 4 matrix that agrees with the number of elements stored in the ExacTrac configuration file matrices.
One last technical issue in rendering is referred to as the w c condition, which states that the homogeneous coordinates of the rendered object must be positive. If not, then the imaging plane and object points are on opposite sides of the camera. A simple way to verify the w c condition is to act on the linac isocenter with M and check the sign of the homogeneous element. The result is to multiply M by the sign of its last entry, in the fourth row and fourth column, in order to satisfy the w c condition.
We now have enough information to work backwards from the ExacTrac renderer matrices to extract the geometry for DRR generation. Consider the 3 × 3 submatrix of M consisting of its first three rows and first three columns. The above analysis implies that the submatrix is a product of projection and rotation transformations. We created a function for RQ decomposition of any real square matrix, where R denotes an upper-triangular projection matrix and Q denotes an orthogonal rotation matrix. The ExacTrac renderer matrices in our system are already denormalized, so that the S NDC matrix may be removed from the equation for M. The relevant detector-scaled projection is which corresponds to R in the RQ decomposition. SID information in the 3 × 4 matrix format is not lost, however. ExacTrac saves the homogeneous coordinate in the third row. Consequently, we need to divide P det by its last entry in the third row and third column. We easily read out the SID as SID = P The fourth column of L 1,2 , besides the homogeneous element, is precisely the fourth column of M. Meanwhile, the 3 × 3 submatrix of L 1,2 (equivalently, of R 1,2 ) is precisely Q in the RQ decomposition. Therefore, the vectors t 1,2 are easily derived by inverting a system of linear equations, an operation that was performed in MATLAB. The vectors must be converted from IEC to HFS coordinates. The direction cosines matrices D 1,2 are given by the transpose of Q. ITK requires a DRR origin in HFS coordinates, which we calculate for each DRR as Our repository includes a modified version of the DRR generator available in ITK. 8 For completeness, we used the exact parameters in our custom-made script in order to compare landmark points of interest. Although our original software is much slower than the ITK functions, we include it in this report as an independent verification of ITK ray tracing and as a tool for understanding DRR generation from first principles. Once again, a linear mapping between CT number and voxel intensity is assumed. We assumed an X-ray beam effective energy of 70 keV, which corresponds to a linear attenuation coefficient of 0.029 cm −1 for water. We also fixed a minimum threshold of 100 HU, which was sufficient for our purpose of replicating DRRs of bony anatomy. Thus, any tissue in a voxel with CT number below the threshold was neglected when calculating radiological paths with line integrals. Other attenuation coefficients and thresholds may be more suitable for certain applications. The ray tracing function provided by ITK, however, only integrates voxels along a ray, above a certain intensity threshold, and does not strictly calculate radiological paths. The 6D registration parameters determined by ExacTrac are optional arguments in the generator.
In order to quantify the matching accuracy between both algorithms and the ground truth, we performed a landmark analysis in MATLAB. Corresponding landmark points were identified for stereoscopic DRRs obtained at slightly different views. Images of an SBRT phantom and a patient spine acquired in our clinic were used in the landmark analysis. For each image, we displayed the known DRRs and our DRRs side by side. One landmark point was selected across the three DRRs before proceeding to the next point. Markers for the F I G U R E 2 Digitally reconstructed radiographs (DRR) generation results for a stereotactic body radiotherapy (SBRT) phantom (top row) and a patient treated at our clinic (bottom row), with values derived from ExacTrac using our methodology. (a) Known DRRs created by ExacTrac. (b) DRRs generated by our customized Insight Segmentation and Registration Toolkit (ITK) script. (c) DRRs generated by our ray tracing algorithm built from scratch points already selected were visible throughout the process to ensure that a distinct set of points were identified. For a given algorithm, we defined landmark displacements as the difference in pixel locations between the landmarks from the generated DRR and the Exac-Trac DRR. Finally, we performed rigid image registrations using the elastix software package 13 for an independent assessment of systematic offsets between our DRRs and treatment-position reference DRRs from the ExacTrac system. Rigid registrations were implemented with a multiresolution Euler transform,which depends on three global parameters, namely vertical and horizontal translations and a rotation about the center of the image.

RESULTS
See Figure 2 for a verification that our modified ITK script and our original MATLAB software reproduce the treatment-position reference DRRs of one SBRT phantom and one patient scanned in our clinic, given the ExacTrac MLinToFlat1 and MLinToFlat2 strings and the analysis described here. Statistics that summarize the landmark displacements, recorded in mm, across reconstruction algorithms are provided in Table 1. Specifically, the mean, standard deviation, median, interquartile range (IQR), minimum, and maximum landmark displacements between each generator and the known DRRs were calculated. Mean displacements are 1.15 mm on average. Parameters for rigid registration to the known DRRs are shown in Table 2. Translations are listed in mm, and the counterclockwise angle is listed in deg about the image center. Our results indicate that rotations are negligible, and absolute value vertical and horizontal translations are 0.18 and 0.35 mm on average, respectively.

DISCUSSION
Perspective projections are best explained as a sequential application of linear operators on ℙ 3 . Given an object in an ambient space with fixed origin, the goal is to find an operator that sends points in the object to an imaging plane along pinhole camera rays. First, we initialize the camera position through a translation matrix and the orientation through a rotation matrix. After such a change of basis, we apply a standard calibration matrix that projects voxels to a fixed orthogonal distance away from the camera. In other words, 3D points within an object are mapped to a 2D imaging plane. A detector scaling operator is required to achieve the desired DRR resolution and upper left corner DRR origin. In general, rendering is restricted to a frustum, so that points outside the frustum are clipped. An NDC conversion maps the frustum to a unit cube. Dividing out by the homogeneous, coordinates yield DRR pixel locations. Offline stereoscopic DRR production is a crucial step in independent image analysis, as the ExacTrac system currently does not enable DRR export for user-defined couch repositioning. For example, deep learning registration methods to detect vertebral body misalignments will require the ability to proficiently create DRR training sets with accurate 6D couch motion. Our provided materials, together with inputted ExacTrac configuration files and treatment logs, allow for an efficient radiograph reconstruction that may be incorporated in further image processing. It may be possible to modify other widely available DRR generators using our methods, including plastimatch, 14 DRRGenerator, 15 and DeepDRR. 16 We found that implementing our methods in ITK was the most straightforward, although in principle one could adapt these generators to include additional parameters necessary for ExacTrac DRR reconstruction. Our DRR algorithms were not designed to handle artifacts such as streaking, which other generators may be able to accommodate. Although it is preferrable to have artifactfree DRRs, artifacts are acceptable when building models for misalignment detection research so long as the presence or absence of artifacts is consistent in both the model training and inference data. A generator such as DeepDRR may be able to provide more realistically simulated radiographs, but our primary task was replicating the geometry for applications to misalignment detections.
Algorithms for offline DRR rendering that successfully reproduce ExacTrac DRRs were derived from first principles. The landmark analysis and rigid registration results demonstrate a high level of accuracy between the DRRs generated by our two algorithms and the known DRRs. All parameters necessary for the DRR generators are stored in the ExacTrac configuration files produced for each treatment session. We illustrate excellent agreement between known DRRs and both our custom-made and ITK-based programs on three example scans. A supplemental repository containing the DRR generator materials is available. Given ExacTrac settings for a particular treatment, we explain the steps to derive projective geometry parameters. Our original software highlights the details on how the parameters are used to produce DRRs, while the customized ITK software yields equivalent results, albeit with sharper images and a faster runtime. Our results are a crucial step in IGRT applications such as vertebral body misalignment detection. Developments in such fields may eventually improve the online performance of ExacTrac to further reduce radiotherapy errors for greater patient outcome.

AC K N OW L E D G M E N T S
We would like to thank our correspondence at Brainlab for their insight on perspective projections and 2D rendering. Our understanding of the configuration data stored by ExacTrac would not have been possible without their guidance and support. We especially thank Dr. Birte Domnik who provided the initial documentation for rendering methods used by ExacTrac,including information on normalized calibrations.

C O N F L I C T O F I N T E R E S T
The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

AU T H O R C O N T R I B U T I O N S
All the authors conceptualized and designed the project and also participated in manuscript drafting and final approval of the manuscript. John Charters developed the mathematical methods, gathered the data, and wrote the computer code.