Cornel Zachiu
email: c.zachiu@umcutrecht.nl

Baudouin Denis De Senneville

Rob H N Tijssen

Alexis N T J Kotte

Antonetta C Houweling

Linda G W Kerkmeijer

Jan J W Lagendijk

Chrit Moonen

Mario Ries

Chrit T W Moonen

Non-rigid CT/CBCT to CBCT registration for online external beam radiotherapy guidance

Keywords: Image-guided radiotherapy, CT/CBCT registration, Online guidance

Introduction

Image-guided external beam radiotherapy (EBRT) provides the possibility of radiation cancer treatment delivery with a high degree of geometric precision and accuracy [START_REF] Gupta | Image-guided radiation therapy: Physician's perspectives[END_REF]. This is facilitated by the capability of the current treatment machines to provide onboard 3D imaging of the target pathology and its surroundings. Due to its cost effectiveness and typically low radiation dose to the patient, cone beam computed tomography (CBCT) has been integrated for therapy guidance purposes, in most state-of-the-art therapeutic linear accelerator (linac) solutions [START_REF] Guckenberger | Image-guided radiotherapy based on kilovoltage cone-beam computed tomography -A review of technology and clinical outcome[END_REF]. During EBRT treatments, a CBCT image is usually acquired at the start of a radiation fraction and compared to a high-resolution computed tomography (CT) image acquired during the planning phase of the treatment [START_REF] Thilmann | Correction of patient positionning errors based on in-line cone beam CTs: Clinical implementation and first experiences[END_REF]. The outcome of this comparison is a set of displacements indicating the anatomical situation at the time of treatment relative to therapy planning. These displacements are currently used to reposition the patient such that the pathology at the time of treatment is approximately at the same location as during planning [START_REF] Thilmann | Correction of patient positionning errors based on in-line cone beam CTs: Clinical implementation and first experiences[END_REF][START_REF] Oldham | Cone-beam-CT guided radiation therapy: A model for on-line application[END_REF][START_REF] Lehmann | Commissioning experience with cone-beam computed tomography for image-guided radiation therapy[END_REF]. Also, in case of high fraction doses or long treatment times, an additional CBCT image is acquired after the radiation fraction delivery, in order to determine whether significant motion occurred during treatment delivery. However, depending on the imaged anatomy and rotation speed of the linac gantry, the acquisition of both CT and CBCT images may be hampered by physiological motion [START_REF] Keall | The management of respiratory motion in radiation oncology report of AAPM Task Group 76[END_REF][START_REF] Marchant | Reduction of motion artefacts in on-board cone beam CT by wapring of projection images[END_REF]. Respiratory motion, in particular, can induce large displacements of the organs situated in the thorax and upper abdomen.

If not taken into consideration, this can lead to motion-induced blurring effects in the acquired images, increasing the uncertainties related to the location of the pathology and organs-at-risk [START_REF] Keall | The management of respiratory motion in radiation oncology report of AAPM Task Group 76[END_REF][START_REF] Li | Motion correction for improved target localization with on-board cone-beam computed tomography[END_REF]. In this sense, a common practice is to continuously acquire X-ray projections over a longer period of time, as the linac gantry rotates around the patient, and reorder them according to the respiratory phase in which they were acquired [START_REF] Sonke | Respiratory correlated cone beam ct[END_REF], Dietrich et al. 2006, Jia et al. 2012). The reordered X-ray projections are then used to reconstruct a 4D (3D + time) image series, providing the configuration of the imaged anatomy over the respiratory cycle. Estimation of respiratory displacements based on a 4D CT/CBCT series could be useful, for example, for providing patient-specific motion models which, in turn, can be used to adapt treatment margins on a case-by-case basis [START_REF] Coolens | A margin model to account for respiration-induced tumous motion and its variability[END_REF][START_REF] Dhou | 3D fluoroscopic image estimation using patient-specific 4DCBCT-based motion models[END_REF][START_REF] Cai | 4D cone beam CT-based dose assessment for SBRT lung cancer treatment[END_REF]. This usually results in a more conformal radiation dose for the target pathology while at the same time improving the sparing of healthy tissues. The displacements between CT and CBCT images acquired during the different phases of EBRT can be estimated using a process called image registration [START_REF] Mani | Survey of medical image registration[END_REF]. The output of the registration process between two images is a transformation that allows mapping one of the images, called the moving image, to the other image, called reference, leading to a spatial consistency between the two images. A different number of degrees-of-freedom can be imposed on the estimated transformation, which is usually chosen depending on the mobility of the imaged anatomical site. A rigid transformation which has a low number of degrees-of-freedom is, for example, more suitable for static anatomies such as the brain or parts of the skeletal system [START_REF] Hill | Medical image registration[END_REF]. For mobile and/or deformable organs, on the other hand, a rigid transformation may be sub-optimal. Thus, the number of degrees-of-freedom of the estimated transformation may have to be increased in order to capture the complex deformation underwent by these organs. Regardless of the nature of the estimated transformation, however, registration of CBCT images is a challenging task, especially for registration methods which rely on image intensity. This is due to the fact that, in order to reduce the radiation dose related to imaging and the required acquisition time, CBCT images are acquired/reconstructed using a low number of X-ray projections and/or low beam intensities. This, in turn, leads to the images containing so-called streaking artifacts while at the same time being affected by a low contrast and a low signal-to-noise ratio (SNR) [START_REF] Ahmad | Optimization of 4D cone-beam CT: Evaluation of streaking artifacts and noise with various simulated gantry rotation speeds[END_REF][START_REF] Schulze | Artefacts in CBCT: A review[END_REF]. Streaking artifacts become especially problematic when CT-CBCT intensity-based registration is of interest, due to the intrinsic intensity inconsistencies they introduce between the CT and the CBCT image [START_REF] Li | Evaluation of deformable image registration for contour propagation between CT and cone-beam CT images in adaptive head and neck radiotherapy[END_REF]. Previous studies address this issue by incorporating intensity correction/matching procedures into existing intensity-based registration algorithms such as demon registration [START_REF] Nithiananthan | Demons deformable registration of CT and cone-beam CT using an iterative intensity matching approach[END_REF][START_REF] Zhen | CT to cone-beam CT deformable registration with simultaneous intensity correction[END_REF][START_REF] Park | Deformable registration of CT and cone-beam CT with local intensity matching[END_REF] or optical flow [START_REF] Ostergaard | Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy[END_REF][START_REF] Li | Evaluation of deformable image registration for contour propagation between CT and cone-beam CT images in adaptive head and neck radiotherapy[END_REF]. The efficiency of such approaches, however, depends on the particularities of the artifacts hampering the CBCT images on a case-by-case basis. More precisely it depends on the type, amount and spatial distribution of the artifacts. Thus, they may provide reliable results only in a limited amount of cases. Even though CT and CBCT images are acquired using the same modality, previous studies have also addressed CT-CBCT registration using methods that are typically used for aligning images acquired with different modalities. Due to their functioning principle, such methods relax the assumptions made on image intensity and align the images based on other criteria such as mutual information (MI) or normalized cross correlation (NCC) [START_REF] Greene | Constrained non-rigid registration for use in image-guided adaptive radiotherapy[END_REF][START_REF] Rubeaux | Evaluation of non-rigid constrained CT/CBCT registration algorithms for delineation propagation in the context of prostate cancer radiotherapy[END_REF][START_REF] Lou | Joint CT/CBCT deformable registration and CBCT enhancement for cancer radiotherapy[END_REF]). However, both MI and NCC are evaluated on a global level. Therefore, while such approaches may perform well for estimating rigid displacements, they may become sub-optimal for elastic deformations. Moreover, such methods typically imply high computational demands and optimization of several input parameters, limiting their use in a clinical setting with the patient still on the interventional table. A different approach to CT-CBCT registration involves aligning the images based on a set of common features/landmarks. It is, for example, current clinical practice that several CT/CBCT-visible fiducial markers are implanted inside the tumor and its proximity [START_REF] Kothary | Percutaneous implantation of fiduacial markers for imaging-guided radiation therapy[END_REF][START_REF] Timmerman | Image-guided and adaptive radiation therapy[END_REF][START_REF] Marchant | Automatic tracking of implanted fiducial markers in cone beam CT projection images[END_REF]. The markers are then either automatically or manually identified and re-aligned in order to estimate a displacement between the images. Besides the implantation procedure involving discomfort for the patient (since in some cases it is performed without anesthesia), this only allows a rigid alignment between the images. Other methods which rely on naturally occurring anatomical markers [START_REF] Paquin | Multiscale registration of planning CT and daily cone beam CT images for adaptive radiation therapy[END_REF][START_REF] Kearney | Automated landmark-guided deformable image registration[END_REF] are sensitive to the number and the discriminative properties of the landmarks. The number of landmarks also limits the degrees-of-freedom of the estimated transformation. In addition, both implanted and natural markers are only surrogates for the tumor position. There is a possibility that the tumor shifts relative to the markers, leading to misalignments between the images.

In the current study, the issue of both CT-CBCT and CBCT-CBCT registration is addressed via the EVolution registration algorithm, recently proposed by Denis de Senneville et al in (Denis de [START_REF] Denis De Senneville | EVolution: an edge-based variational method for non-rigid multi-modal image registration[END_REF]. In simplistic terms, the EVolution method estimates a deformation between two images, by locally matching similar contrast patterns. Thus, contrary to intensity-based methods, the algorithm is by construction robust to local and/or global intensity changes from one image to the other. Also, it does not require the implantation, detection and/or identification of a set of features/markers. At the same time it is capable of providing an elastic transformation between the registered images on a voxel-by-voxel basis. The method requires a small number of input parameters and implies a fast numerical scheme rendering it as an attractive solution for online therapy guidance. The contribution of the current work is fourfold: (i) A validation procedure is proposed, which allows quantifying the impact of the artifacts and low SNR of the CBCT images on the motion estimates.

(ii) The performance of the EVolution algorithm is evaluated qualitatively and quantitatively for CT to CBCT registration.

(iii) The algorithm's performance for CBCT-based 3D respiratory motion estimation was investigated.

(iv) An implementation with low computational demands is performed for the EVolution method, making it suitable for use in a clinical setting for online therapy guidance. et al. 2016) and it provides the displacement between two images as the optimizer of the following functional:

Method description

E(T ) = Ω e -C(T ) + α 2 ∇u 2 2 + ∇v 2 2 + ∇w 2 2 d r (1)
with

C(T ) = Γ ∇ I (T ( r)) • ∇ J ( r) d r Γ ∇ I (T ( r)) 2 ∇ J ( r) 2 d r (2)
where T = (u, v, w) is the 3D displacement, I and J are the two images to be aligned, Ω is the image domain, Γ is a cubic local neighborhood around the voxel being processed, r is a spatial location, ∇ is the 3D gradient operator, • 2 is the Euclidean norm and α is an algorithm parameter whose purpose will be detailed shortly. Notice that the functional optimized by the EVolution algorithm (Eq. 1) is defined by two terms. The first term of the integral, also known as the data fidelity term (or simply data term), is a measure of the similarity between the images. The rationale behind using this term as a similarity metric becomes apparent if C(T ) is re-written under the following form:

C(T ) = Γ w T ( r) |cos(∆θ T ( r))| d r Γ w T ( r)d r (3)
with w T ( r) and ∆θ T ( r) given by:

w T ( r) = M I (T ( r))M J ( r) ∆θ T ( r) = θ I (T ( r)) -θ J ( r) (4)
where M I and M J are the magnitudes of the spatial gradients of the two images and ∆θ is the phase difference between the gradient orientations. Therefore, the term w T ( r) favors the alignment of edges/gradients that are present in both images, while |cos(∆θ)| provides a measure for the angle between the gradient orientations, favoring parallel and anti-parallel gradients and implicitly similar local contrast patterns. The optimization of the data term alone is generally an ill-posed problem with an infinity of solutions. This was addressed by adding a regularization term to the functional in Eq. 1, which constrains the estimated deformation to be spatially smooth/differentiable. The amount of smoothness is controlled by the parameter α: the larger the value of α, the smoother the resulting deformation will be.

2.1.2. Numerical scheme In the current work, the functional in Eq. 1 was optimized using an approach similar to the one proposed in the seminal paper of Denis de Senneville et al. In effect, its minimizer was found by solving the associated Euler -Lagrange equations. This ultimately implied finding the solution of a non-linear system. In order to reduce the non-linearity of the functional, a coarse-to-fine strategy was employed, which iterated the algorithm from a 32-fold downsampled version of the images step-by-step to the original resolution. The motion field estimated at each step was used as an initialization for the next higher resolution level. At each resolution level, the algorithm was considered to have converged when the average difference between the displacements estimated at the current and the previous iteration was smaller than 10 -3 voxels. For a more detailed and technical description of the numerical scheme, please consult (Denis de Senneville et al. 2016).

Experimental setup

In order to test the capability of the EVolution algorithm for CBCT-based motion estimation, several experiments were conducted:

(i) The performance of the EVolution algorithm was initially evaluated for CT to CT registration and put in perspective with respect to a state-of-the-art method (namely the optical flow algorithm [START_REF] Zachiu | An improved optical flow tracking technique for real-time MR-guided beam therapies in moving organs[END_REF]). This evaluation was of interest due to the fact that CT to CT registration played an important role for quantifying the accuracy and precision of the EVolution method for both CT-CBCT and CBCT-CBCT registration.

(ii) A validation procedure was developed, based on CBCT images synthesized from high resolution CT images. Accuracy and precision of the method was then evaluated for CT -synthetic CBCT registration, with the synthesized CBCT containing an increasing amount of streaking artifacts.

(iii) A qualitative evaluation of the algorithm was performed for registering CT to CBCT images, each acquired with their respective imaging hardware.

(iv) Finally, the algorithm's performance for 3D respiratory motion estimation based on synthesized CBCT images was investigated.

The following sections will provide details on each of the 4 experiments.

2.2.1. Validation of the EVolution algorithm for CT to CT registration As it will be detailed in the following sections, CT to CT registration plays an important role, in the current work, for quantifying the accuracy and precision of the EVolution algorithm for CT to CBCT and CBCT to CBCT registration. For this reason, a preliminary experiment was conducted on CT data acquired on 5 head-and-neck cancer patients. Each data set contained a pair of CT images acquired at different time instants (more than 1 week apart). On each of the CT images, an experienced radiation therapist delineated both pathological areas and organs-at-risk (OAR). The 5 pairs of CT images were then registered using both EVolution and the optical flow algorithm [START_REF] Zachiu | An improved optical flow tracking technique for real-time MR-guided beam therapies in moving organs[END_REF] (used for comparison). For both registration methods, the Dice Similarity Coefficient (DSC) [START_REF] Li | Comprehensive evaluation of ten deformable image registration algorithms for controu propagation between ct and cone-beam ct images in adaptive hean & neck radiotherapy[END_REF], normalized cross correlation (NCC) [START_REF] Luo | A fast normalized cross-correlation calculation method for motion estimation[END_REF], structural similarity index (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] were evaluated after registration for several anatomical sites. In addition, the voxel-wise error in flow endpoint (FEP) between the EVolution and the optical flow motion fields was evaluated [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]:

F EP ( r) = u E Vo( r) -u OF ( r) 2 (5)
where u EVo = (u EV o , v EV o , w EV o ) and u OF = (u OF , v OF , w OF ) are the 3D motion fields estimated by the EVolution and the optical flow algorithms, respectively, • 2 is the Euclidean norm and r is the voxel position. For computational purposes, the images and organ delineations were projected on a 256 × 256 × 256 lattice with a voxel size of 1.5 × 1.5 × 1.5 mm 3 . The axes of the lattice were aligned with the axes of the DICOM patient coordinate system (NEMA PS3 / ISO 12052 2017).

Validation of the EVolution algorithm using synthetic CBCT images

Gold standard selection In order to evaluate the precision and accuracy of the EVolution algorithm for CT to CBCT and CBCT to CBCT registration, the selection/establishment of a ground truth deformation was necessary. In general, obtaining an in-vivo gold standard is a challenging task, especially for deformable/elastic anatomies. In the current work, a particular validation strategy was developed in this sense. In a first step, two high resolution CT images acquired at different time points are registered to one another using the EVolution algorithm, with the resulting displacements being established as gold standard. One or both of the images are then used to synthesize CBCT volumes (the synthesis procedure will be described in the next section) and then the images are registered to one another a second time. The resulting displacements are then compared to the gold standard. For CT to CT registration, the EVolution algorithm was already put into perspective with respect prior art (see section 2.2.1). Therefore, this process will provide the individual impact the undersampling of the sinogram has on CT-CBCT and CBCT-CBCT registrations performed using EVolution.

Synthesis of CBCT images from high-resolution CT volumes In order to synthesize a CBCT from a high resolution CT image, the freely available TIGRE CBCT reconstruction toolbox was employed [START_REF] Biguri | TIGRE: A MATLAB-GPU toolbox for CBCT image reconstruction[END_REF][START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. The software allows a basic simulation of a CBCT acquisition and reconstruction, using any image as input, while at the same time providing access to the associated acquisition and reconstruction parameters. In particular, this allows control over the number of X-Ray projections used to reconstruct the CBCT volume. Therefore, the quality of the motion estimates as a function of the number of projections of can be investigated.

Performance evaluation of the EVolution algorithm for CT to synthetic CBCT registration The capability of the EVolution algorithm for registering CT to synthesized CBCT images was evaluated on data acquired on a head-and-neck cancer patient. The data initially consisted in a pair of CT images acquired on different occasions. In a first step towards the evaluation, the dependency of the estimation errors on the amount of the streaking artifacts present in the synthetic CBCT image was investigated. For this purpose, one of the CTs was used to synthesize a CBCT image with an increasing amount of streaking artifacts. The CBCT image was reconstructed using 360, 120, 72, 50 and 36 projections, which corresponds to an acquisition with an angular increment of 1 • , 3 for details on the gold standard). The comparison between the gold standard and the estimated motion field was performed in terms of the spatial distribution of the voxelwise error in flow endpoint (FEP):

F EP ( r) = u( r) -u gold ( r) 2 (6)
where u = (u, v, w) is the estimated 3D displacement, u gold = (u gold , v gold , w gold ) is the gold standard 3D displacement. Following this initial part of the experiment, the performance of the algorithm was investigated in more detail for a CBCT image synthesized using 72 projections. This case was of particular interest, since a reconstruction using 72 projections provides an image quality towards the lower end of what a typical CBCT system would provide. For computational purposes the images were projected, prior to registration, on a 256 × 256 × 256 lattice with a voxel size of 1.5 × 1.5 × 1.5 mm 3 . The axes of the lattice were aligned with the axes of the DICOM patient coordinate system.

CT to clinical CBCT image registration

This experiment aimed to evaluate the algorithm's capabilities for registering CT to clinical CBCT images. In effect, the algorithm was employed for the registration of a high resolution planning CT image to one of the daily CBCT images. This was carried-out for two pairs of CT -CBCT images: one acquired on a kidney cancer patient and the other on a lung cancer patient. The performance of the algorithm was assessed both by visual inspection and quantitative evaluation. The latter was achieved by manual annotation of 10 landmarks in both the kidney and the lung data sets. The resulting manually determined displacements were then compared in terms of the FEP to the displacements estimated by the EVolution algorithm. As an additional registration error metric, the normalized cross correlation (NCC) between the images before and after registration was also calculated.

Before registration, all the images were projected on a 256 × 256 × 256 lattice with a voxel size of 0.9 × 0.9 × 0.9 mm 3 . This was carried-out not only for computational purposes, but also, more important, due to the fact that the original CT and the CBCT images were reconstructed with different voxel sizes. Moreover, for both patients, part of the CT images was masked-out since their field-of-view was considerably larger than the one of the CBCTs, which increases the risk of the algorithm converging towards a local minimum.

2.2.4. Synthetic CBCT-based respiratory motion estimation using EVolution A 4D CT dataset was acquired on a renal cell carcinoma patient. The images sample the anatomy of the upper abdomen at 10 phases of the patient's respiratory cycle, between 0% and 90% of the cycle with a 10% increment. The 4D series was registered using the EVolution algorithm to a reference position provided by the CT image acquired at the 0% phase of the respiratory cycle. The resulting displacements were then stored and used as a gold standard. A 4D CBCT series was then synthesized based on the original CT series, and registered to the CBCT image corresponding to the 0% respiratory cycle phase. The resulting displacements were then compared to the aforementioned gold standard, in terms of the temporally averaged FEP. Additional validation was performed by analyzing the CBCT-estimated 3D trajectory of two landmarks over the respiratory cycle. The first of the landmarks was one of the fiducial markers implanted in one of the kidneys, while the second marker corresponded to the centroid of the contralateral kidney. For the implanted marker, the ground truth was considered to be the trajectory provided by the EVolution algorithm on the original 4D CT series. However, for the centroid of the contralateral kidney, the trajectory estimated on the synthetic 4D CBCT series was compared to a manually determined trajectory. The estimated and the gold standard trajectories were again compared in terms of the FEP for the selected landmarks. The 4D CBCT images were synthesized using 120 and 72 projections respectively, which produced images of average to low quality compared to what most clinical systems provide. Three scenarios were investigated: 1) All CBCT images were synthesized using 120 projections, 2) The reference image was reconstructed using 120 projections while the rest of the images were synthesized using 72 projections and 3) All CBCT images were reconstructed using 72 projections. For computational purposes, the original 4D CT images were downsampled from 512 × 512 × 123 to a size of 256 × 256 × 123 with a 2 × 2 × 3 mm 3 voxel size.

Hardware, implementation and algorithm configuration

The EVolution algorithm was implemented on a graphical processing unit (GPU) using the compute unified device architecture (CUDA) platform. The implementation was then executed on an nVidia Tesla K20 graphics card. EVolution requires as input two parameters, namely α and Γ (see Eq. 1 and 2). A procedure was developed in the scope of the current work in order to calibrate the two parameters. Three landmarks were manually tracked in the 4D CT series described in section 2.2.4. These landmarks consisted in two of the fiducial makers implanted at different locations in the ipsilateral kidney and the centroid of the contralateral kidney. The 4D series was registered to the CT image acquired at the 0% phase of the respiratory cycle, while varying the value of α in the 0.05 -1.0 range, with a 0.05 increment and the size of Γ in the 3 × 3 × 3 -21 × 21 × 21 range with an increment of 2 × 2 × 2. The combination of α and Γ which provided the smallest average FEP between the estimated and the manually determined displacements for the three landmarks was chosen for use.

A similar procedure was carried-out in order to calibrate α and Γ for CT -CBCT registration, with the difference that the reference image was replaced with a synthetic CBCT volume, generated from the CT image acquired at the 0% phase of the respiratory cycle. The CBCT image was synthesized using 72 projections.

For the optical flow algorithm, used in the experiment described by section 2.2.1, an implementation similar to the one described in [START_REF] Zachiu | An improved optical flow tracking technique for real-time MR-guided beam therapies in moving organs[END_REF] was carried-out.

The optical flow algorithm provides motion estimates as the minimizers of the following functional:

E OF ((u)) = r∈Ω (I( r) -J( r + u( r))) 2 + β 2 ∇u( r) 2 2 (7)
where I and J are the reference and the moving image, (u) = (u, v, w) is the 3D set of displacements, r is a spatial location, Ω is the image domain, ∇ is the gradient operator,

• 2 is the Euclidean norm and β is a parameter linking the two terms of the functional. Details related to the numerical scheme used to solve Eq. 7 can be found in [START_REF] Zachiu | An improved optical flow tracking technique for real-time MR-guided beam therapies in moving organs[END_REF]. Similar to α and Γ, the input parameter β was calibrated via an exhaustive search for the value which provided the smallest FEP between the estimated and the manually determined displacements of the same three landmarks. During the search, the value of β was varied in the 0.01 -1.0 interval, with an increment of 0.05 (except for the first incrementation which was 0.04).

Results

Section 3.1 provides an evaluation of the EVolution algorithm's capability for registering high resolution CT images. In addition, the associated results are put in perspective with respect to the optical flow algorithm. Sections 3.2 and 3.3 showcase the algorithm's performance for registering CT to CBCT images for both synthetic and clinical CBCT images. The potential of the method to perform CBCT-based 3D respiratory motion estimation is showcased in section 3.4, followed by a report on the computational requirements of the EVolution method in section 3.6.

Assessment of the EVolution algorithm for CT to CT registration

Table 1 reports the DSC, NCC and SSIM for several volumes-of-interest, after registering pairs of CT images acquired at different time instants on 5 head-and-neck cancer patients. Registration was performed using both the EVolution and the optical flow algorithm. Note that, with very few exceptions, the DSC for both methods remains above 0.7 for all anatomies under consideration, with a consistent value of ∼0.98 for the body contour itself. Moreover, it can be observed that in the majority of cases, the two registration algorithms provide comparable DSC values. Concerning the NCC and SSIM global criteria, the two methods again provide similar results. Table 2 reports the FEP between the motion fields provided by the EVolution and the optical flow algorithms. Evaluation was performed for the same CT data analyzed in Table 1. In the majority of cases, the mean and the standard deviation of the FEP between the two methods remains sub-millimeter. Note that the higher FEP values in

Algorithm performance for CT to synthesized CBCT registration

The accuracy and precision of the EVolution algorithm following the alignment of a CT to synthesized CBCT images was evaluated for a head-and-neck cancer patient (more precisely Patient #1 from Table 1). Methodological details related to the evaluation can be found in section 2.2.2. Fig. 2(a) -2(e) illustrate a transverse slice from the synthetic CBCT images. More precisely, the images were synthesized using 360, 120, 72, 50 and 36 projections. From a visual analysis it can be observed that as the number of projections decreases, the streaking artifacts in the synthesized CBCT images become stronger. Fig. 2(f) -2(j) on the other hand, showcase the spatial distribution of the FEP associated to each of the images from Fig. 2(a) -2(e). Notice that, especially in homogeneous areas, the alignment errors have a tendency to increase with the amount of artifacts present in the images, coming close to 5 mm when 36 projections are used for the reconstruction of the CBCT image. The registration errors were analyzed in more detail for a CBCT image reconstructed using 72 projections (see section 2.2.2 for details). Fig. 3(a) -3(c) displays a coronal, a sagittal and a transverse slice selected from this CBCT image. The red overlay on the three images corresponds to the location of the clinical tumor volume (CTV). Fig. 3(d) -3(f) illustrate the spatial distribution of the FEP prior to registration. This was evaluated by setting u( r) in Eq. 6 to zero. Since u gold ( r) is assumed to be the true motion, by setting u( r) equal to zero the FEP provides the misalignments that would occur in the absence of a registration procedure. It can be observed that, locally, the alignment errors come close or even exceed 25 mm. The spatial distribution of the FEP after registration is illustrated in Fig. 3(g) -3(i). With a few exceptions, the alignment errors rest well beneath 1.5 mm, which corresponds to the image voxel size.

A statistical analysis of the FEP for the pathological tissues and several organsat-risk was also performed. This is reported, before and after registration, in the second and third columns of Table 3. The error statistics are provided under the format mean±standard deviation of a set that includes the pooled FEP in all the voxels contained by the analyzed anatomical structures. On average, after registration, the alignment errors remain sub-voxel. The same observation can be made for the precision of the estimated motion, indicated by the FEP standard deviation. The DSC values before and after registration are reported in columns 4 and 5 of Table 3. Improvements can be observed after registration, for all VOIs, with DSC values ranging from 0.7 for the smaller structures such as the CTV, up to 0.98 for the body itself. An evaluation of the NCC before and after image registration also showcases an increase from 0.89 to 0.96. the CBCT image, a slice from the CT image, checkerboard overlap between the images before registration and checkerboard overlap between the images after registration. An improvement in image alignment can be observed, especially for the implanted fiducial markers (whose coherence is restored after registration), for the skeletal structures and at organ boundaries. As specified in section 2.2.3, 10 landmarks were manually identified in both the CBCT and the CT image. Comparing the manually determined displacements for the landmarks to the ones provided by the EVolution algorithm lead to the observation that the FEP in the 10 landmarks was reduced from 9.57 ± 1.61 mm in the absence of registration, to 1.18 ± 0.3 mm after registration. The alignment errors in the absence of registration were computed by setting u in Eq. 6 to 0. In addition, after registration, the NCC showcased an increase from 0.86 to 0.94. A set of images following the same arrangement as the ones in Fig. 4 is illustrated in Fig. 5, however, the kidney case was replaced with images acquired on a lung cancer patient. Similar to the kidney case, a notable improvement in image alignment can be observed after registration, especially in the lung parenchyma, the skeletal structures and the organ boundaries. Manual tracking of 10 landmarks revealed a decrease of the alignment errors from 6.58 ± 1.54 mm in the absence of registration, to 0.94 ± 0.3 mm with registration. After registration, the NCC also increased from 0.81 to 0.96.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)
(a) (b) (c) (d) (e) (f) (g) (h) (i)

Synthetic CBCT-based 3D respiratory motion estimation

As detailed in section 2.2.4, for this particular experiment, CBCT images were synthesized using as a base a 4D CT series covering one respiratory cycle of a renal cell carcinoma patient. Fig. 6(a) -6(c) illustrate a coronal, a sagittal and a transverse slice, selected from the CBCT image used as reference during the registration process. Fig. 6(d) -6(f) and Fig. 6(g) -6(i) display for the slices illustrated in Fig. 6 the temporally averaged FEP, before and after registration. It can be observed that errors which exceed 5 mm in moving organs such as the liver and kidneys, are reduced, in most areas, to values smaller than ∼1.5 -2 mm. Note that the results illustrated in Fig. 6 were obtained from a 4D CBCT series synthesized using 120 projections. The performance of the EVolution algorithm for CBCT-based respiratory motion estimation was further evaluated by analyzing the estimated 3D trajectory of two landmarks, over the respiratory cycle. One of the landmarks was a fiducial gold marker (indicated by the red arrow in Fig. 6), while the second landmark was the centroid of the contralateral kidney (delineated by the blue dashed line in Fig. 6). The gold standard trajectory for the fiducial marker was provided by an estimation on the original 4D CT image series, while for the contralateral kidney centroid, manual tracking was performed. Fig. 7 illustrates the trajectory of the two landmarks in the head -foot, anterior -posterior and the left -right directions. It can be observed that, with very few exceptions, the distance between the estimated curves and the gold standard remain sub-millimeter. Note that the CBCT images were synthesized from the original 4D CT series using 120 projections. Table 4 reports the errors between the trajectories estimated for the landmarks and their corresponding gold standard trajectory, for three scenarios: 1) Both the reference and the moving CBCT images were synthesized using 120 projections (second row); 2) The reference image was synthesized using 120 projections, while the moving images were synthesized using 72 projections (third row); 3) Both the reference and the moving images were synthesized using 72 projections (fourth row). The first row in the table reports the FEP for the two landmarks, without registering the 4D synthetic CBCT series. The FEP is reported in the mean ± standard deviation format, with the statistics computed on the pooled FEPs from all the sampled phases of the respiratory cycle. Notice that for the two scenarios in which the reference CBCT image is reconstructed using 120 projections, the errors remain approximately in the sub-millimeter range, regardless of the amount of the streaking artifacts in the moving images. However, when both the reference and the moving images were synthesized using 72 projections, the quality of the motion estimates in the contralateral kidney is affected. 4: FEP between the trajectories estimated for the implanted marker and the contralateral kidney, for different amounts of streaking artifacts added to the synthesized 4D CBCT images. The first row indicates the FEP for the two landmarks, in the absence of a registration scheme. The subsequent rows report the errors between the estimated and the gold standard trajectories for the cases in which both the reference and the moving images were synthesized using 120 projections (second row), the reference image was synthesized using 120 projections and the moving images were synthesized using 72 projections (third row) and when both the reference and the moving images were reconstructed using 72 projections (fourth row).

(a) -6(c) (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
M (a) M (b) M (c) (d) (e) (f) (g) (h) (i)

Input parameter calibration for the EVolution and optical flow algorithm

The EVolution algorithm requires α and Γ as input parameters (see Eq. 1 and 2), while the optical flow algorithm requires the calibration of the regularization parameter β (see Eq. 7). Details related to the manner in which these input parameters were optimized can be found in section 2.3. The dependency between estimation accuracy and the parameters α and Γ for both CT to CT and CT to synthetic CBCT registration via the EVolution algorithm is displayed in Fig. 8. The highest accuracy of 0.63 mm and 0.67 mm, respectively, is achieved for the combination (α, Γ) = (0.6, 11 × 11 × 11) for both calibrations. It can also be observed that the accuracy remains close to its optimal value for a rather wide range of values of α and Γ. A rapid decrease in algorithm accuracy can be noticed as the value of α drops under 0.25 -0.3. Fig. 9 illustrates, in the scope of this study, the relationship between the accuracy of the optical flow algorithm and the regularization parameter β. It can be observed that the best accuracy is ∼1.2 mm, which is achieved for a value of β equal to 0.05 (value which was employed in the present work).

Computational performance of the algorithm

Following the benchmarking of the computational requirements of the EVolution algorithm it was found that for images of size 256×256×256, the algorithm converges on average in approximately 60 s.

Discussion

Several concepts in image guided radiotherapy, such as day-to-day positioning compensation [START_REF] Thilmann | Correction of patient positionning errors based on in-line cone beam CTs: Clinical implementation and first experiences[END_REF][START_REF] Oldham | Cone-beam-CT guided radiation therapy: A model for on-line application[END_REF][START_REF] Lehmann | Commissioning experience with cone-beam computed tomography for image-guided radiation therapy[END_REF], "virtual couch shifts" [START_REF] Bol | Virtual couch shift (VCS): Accounting for patient translation and rotation by online IMRT re-optimization[END_REF], dose accumulation [START_REF] Richter | Investigation of the usability of conebeam CT data sets for dose calculation[END_REF][START_REF] Sykes | Dosimetry of CBCT: methods, doses and clinical consequences[END_REF] or delineation propagation [START_REF] Rubeaux | Evaluation of non-rigid constrained CT/CBCT registration algorithms for delineation propagation in the context of prostate cancer radiotherapy[END_REF]) often rely on establishing a spatial coherence, using image registration, between a planning high-resolution CT and a daily CBCT image and/or between CBCT images acquired at different time instants. In order to limit the radiation dose associated to imaging, which is of particular importance for pediatrics [START_REF] Olch | Pediatric radiotherapy: Planning and treatment[END_REF], the CBCT volumes are frequently reconstructed using a low number of X-Ray projections and/or a low amount of beam intensity. This leads to low SNR images altered by streaking artifacts, which is problematic for state-of-theart registration algorithms relying on gray-level intensity conservation. Moreover, the registration process is also often hampered by the low soft-tissue contrast, which is intrinsic to the CT and CBCT images. The current work demonstrates, that despite such constraints on image quality, the recently proposed EVolution algorithm (Denis de Senneville et al. 2016) is fully capable of providing reliable estimates of the deformation between the images, with a high degree of precision and accuracy. In an initial step, however, prior to investigating the performance of the EVolution algorithm for CT-CBCT and CBCT-CBCT registration, the current study evaluates its capabilities for re-aligning CT images, since CT to CT registration played a key role during the validation procedures employed in the current work (see, for example, section 2.2.2). This was achieved by evaluating the post-registration DSC, NCC and SSIM provided by the EVolution algorithm, for several anatomical structures on CT data acquired on 5 head-and-neck cancer patients. During this experiment, the use of the DSC, NCC and SSIM as registration quality evaluation metrics is argued by their previous use in several independent studies focusing on CT and CBCT image registration [START_REF] Kumarasiri | Deformable image registration vased automatic CT-to-CT contour propagation for head and neck adaptive radiotherapy in the routine clinical setting[END_REF][START_REF] Rigaud | Evaluation of deformable image registration methods for dose monitoring in the head and neck radiotherapy[END_REF][START_REF] Li | Comprehensive evaluation of ten deformable image registration algorithms for controu propagation between ct and cone-beam ct images in adaptive hean & neck radiotherapy[END_REF][START_REF] Park | Deformable registration of CT and cone-beam CT with local intensity matching[END_REF], providing a means to relate the obtained results to prior art. In the current work, the obtained DSC, NCC and SSIM values were also put in perspective with the ones provided by a pre-existing optical flow algorithm. The latter was selected as a point-of-comparison due to its prior evaluation and validation for CT to CT registration in the study of Ostergaard et al [START_REF] Ostergaard | Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy[END_REF]. Therein it was demonstrated that the optical flow algorithm can achieve an accuracy of ∼1.1 mm, value which was concluded to be acceptable for registering CT volumes. Such results are also in good correspondence with the ones obtained in the current work (see section 3.5). In the scope of this study, both the EVolution and the optical flow algorithm provided similar DSC, NCC and SSIM values. The DSC in particular ranged from ∼0.7 for smaller structures up to 0.98 for the body contour. These values are overall comparable with the ones provided by other state-of-the-art methods [START_REF] Kumarasiri | Deformable image registration vased automatic CT-to-CT contour propagation for head and neck adaptive radiotherapy in the routine clinical setting[END_REF]. Also, an evaluation of the FEP between the motion fields provided by the EVolution and the optical flow applied on the same CT datasets, resulted on average in sub-millimeter values (see Table 2). There are, however, a few isolated cases with rather large discrepancies between the DSC values provided by the EVolution and the optical flow algorithm (see, for example, left and right parotid in patients #4 and #5 from Table 1). These discrepancies are also in good correspondence with large FEP values in Table 2. This can be explained by the different fundamental principles on which the two algorithms rely for performing motion estimation. While EVolution assumes the conservation of contrast between the reference and the moving image, the optical flow algorithm assumes voxel intensity conservation. Therefore, depending on which of the two assumptions is locally violated, one or the other algorithm will perform poorer in a particular area. This may also explain the rest of the smaller differences between the DSC values provided by the two compared algorithms. It is also worth noting that, in the current work, validation of the EVolution algorithm for CT to CT registration was performed solely for confirmation purposes, since this aspect was already addressed in the seminal work of Denis de Senneville et al (Denis de [START_REF] Denis De Senneville | EVolution: an edge-based variational method for non-rigid multi-modal image registration[END_REF]. Over the course of the presented work, two quality criteria were employed for evaluating the outcome of image registration: the DSC and the FEP. The reason for which one was chosen over the other depended on the gold standard which was available. For CT to CT registration on the head-and-neck cancer patient data, the only available ground truth were the structure delineations made on the reference and the moving image. This does not allow the calculation of the FEP, since gold standard displacements were not available for these data sets. Therefore, the DSC was chosen as a quality evaluation metric in this case. The downside of the DSC is that it only provides a qualitative measure of registration performance. On the other hand, for the experiments where displacements or point-to-point correspondences were available as a gold standard, the FEP was chosen for use, since it provided a quantitative measure of registration accuracy and precision. Nevertheless, it should be noted that, in case of a severely low number of gold standard displacements or point-to-point correspondences between the reference and the moving image, the DSC may be a better choice (since it usually encompasses a larger number of voxels). In the scope of this paper, the performance of the EVolution algorithm for CT to CBCT registration was evaluated using three gold standards: the motion fields provided by the registration of high resolution CT images, the DSC for several anatomical structures, including pathological tissues and OARs and manually annotated landmarks. Qualitative visual criteria were also employed for the evaluation, in particular for the study described in section 2.2.3. Fig. 2 showcases the impact on the motion estimates of an increasingly coarse sampling of the sinogram for the reconstruction of the analyzed CBCT image. Note that even for CBCT images reconstructed using only 72 projections, leading to an image quality at the lower end of what a typical clinical system would provide, the additional registration errors remain predominantly sub-millimeter. Even for less than 72 projections, only a few error peaks start occurring in the more homogeneous areas of the image. The body contour and skeletal muscle are still correctly re-aligned. These observations were further confirmed during the more in-depth analysis conducted on one of the headand-neck cancer patients. Fig. 3 showcases the spatial distribution of the alignment errors in the aforementioned patient, following the registration of a CT to a CBCT image reconstructed using 72 projections. The errors are illustrated in three planes intersecting approximately midway into the CTV. Noteworthy is that the alignment errors remain mostly sub-millimeter, including for the nasal and air cavities, which have proven problematic in previous related studies [START_REF] Paquin | Multiscale registration of planning CT and daily cone beam CT images for adaptive radiation therapy[END_REF][START_REF] Park | Deformable registration of CT and cone-beam CT with local intensity matching[END_REF][START_REF] Zhen | CT to cone-beam CT deformable registration with simultaneous intensity correction[END_REF]). The only locations in which the errors of the EVolution method have a tendency to increase is in areas devoid of contrast, such as deep into the brain and isolated into the shoulder muscles. There are basically two factors which contribute to this: For one, the lack of intrinsic contrast in the CTs in these areas shifts the burden of the registration process towards the regularization term. Second, since the CBCT artifacts in such areas are dominant compared to ''true'' anatomical contrast, they counteract the effect of the regularization, which as a net-effect degrades the quality of the registration in these domains. Nevertheless, even in these regions, with few exceptions, the errors remain close to the voxel size (∼1.5 mm). The statistical distribution of the alignment errors for the pathology and several of the organs-atrisk of the same patient is reported in Table 3. Following registration, sub-millimeter precision and accuracy was achieved for all analyzed anatomical structures. Also, by comparing the post-registration DSC values reported in Table 3 to the ones in Table 1 (Patient #1) it can be observed that, despite the additional artifacts and low SNR altering the CBCT image, the DSC remains approximately the same (≥0.7 for smaller anatomies up to 0.98 for the body contour) as if a CT acquisition would have been used instead. The good performance of the EVolution algorithm for CT-CBCT registration is further confirmed by the kidney and lung cases illustrated in Fig. 4 and 5, respectively. Following registration, in both cases, a high degree of spatial continuity from one square to the next in the composite images can be observed for the visible features. The most obvious are the implanted fiducials in the kidney cancer patient, the parenchyma of the lung cancer patient and the organ boundaries and skeletal structures in both patients. Comparing the manually determined displacements for 10 landmarks in each of the two data sets to the ones provided by the EVolution algorithm revealed an accuracy and precision of the algorithm comparable to the voxel size. In addition, both cases have demonstrated after registration an increase of the NCC to values which are in good correspondence with the state-of-the-art [START_REF] Park | Deformable registration of CT and cone-beam CT with local intensity matching[END_REF]. Concerning its capability for CBCT-based 3D respiratory motion estimation (i.e. for CBCT to CBCT image registration), the EVolution algorithm was validated against two gold standards: the motion fields resulting from the alignment of the CT images used to synthesize the 4D CBCT series and the trajectory of two landmarks over a respiratory cycle. Fig. 6 displays the spatial distribution of the temporally averaged alignment errors for the 4D CBCT series. In particular for the kidneys, which were the organs of interest in the illustrated data, it can be observed that after the alignment of the 4D series, the errors remain mostly sub-millimeter. On a global scale, except for a few isolated areas, the post-registration misalignments rest under the in-plane voxel size (2 mm). The good capabilities for respiratory motion estimation of the EVolution algorithm are further emphasized in Fig. 7, where it can be observed that the estimated 3D trajectory of the selected landmarks closely follows the gold standard trajectory. A more quantitative evaluation of the errors between the estimated and the gold standard trajectories is provided in Table 4. In the table it is shown that, as long as the sinogram sampling of both CBCT images (reference and moving image) resides above a particular threshold, sub-millimeter precision and accuracy can be achieved. Also, asymmetric registrations with a high-quality reference image and a "faster" CBCT image of lower quality maintains an accuracy and precision below one millimeter. It is only when both the reference and the moving images of the 4D CBCT series are reconstructed using 72 projections that the algorithm performance starts to degrade below the native image resolution. An interesting aspect is that, as the quality of the synthetic CBCT images degrades, only the registration errors for the centroid of the contralateral kidney have a tendency to increase (see Table 4). This can be explained by the presence of the implanted fiducial markers in the ipsilateral kidney, which appear as hyper-intense features in the images. Therefore, despite an increase in the magnitude of the streaking artifacts, the fiducial markers provide sufficient contrast in the synthetic CBCT images, facilitating and accurate and precise tracking. An issue which was generally overlooked by previous related studies is that the CT and the CBCT images are often acquired using different geometries (e.g. FOV, image size, voxel size, etc.). If not taken into account, this can frequently lead to unreliable motion estimates. In the current study this was addressed by regridding the images into a common coordinate system prior to registration. The dimensions and voxel size of the new lattice were chosen such that the entire FOV of the original images is also covered by the regridded images. Obviously, the two lattice parameters (i.e. image and voxel size) can also be optimized with respect to computational performance since, in general, the computational time required by the EVolution algorithm depends on the size of the images. In addition, for the kidney and lung cancer patients (see Fig. 4 and5) the FOV of the CBCT images was considerably smaller than that of the CT images, which is often the case in clinical practice. If not taken into consideration, this may impair the convergence of the algorithm to the global minimum. In the scope of this paper, this was addressed by imposing border conditions on the CBCT image such that all voxels that are outside or on the border of the FOV are excluded from the registration process. One of the traits that makes EVolution an attractive option for motion estimation in a clinical setting is the low number of input parameters it requires: the regularization parameter α and the patch size Γ (see Eq. 1 and 2). A noteworthy fact is that, once Γ and α were optimized, their value was maintained constant during all experiments, with a consistent good quality of the motion estimates. Moreover, both the calibration for CT to CT registration and for CT to synthetic CBCT registration provided the same optimal values for α and Γ (see Fig. 8). So it is fair to say that for the problem at hand, no inter-individual adjustments of the parameters were necessary. Therefore, we hypothesize, that once the acquisition scheme of the CBCT images is fixed, the values of α and Γ could be optimized prior to the intervention, making algorithm configuration completely transparent to the clinician delivering the therapy. One elegant way to optimize the value of the input parameters prior to the intervention is, for example, to register two CT images acquired on a patient and establish the resulting motion field as gold standard. One or both of the images (depending on whether CT-CBCT or CBCT-CBCT registration is of interest) could then be used to synthesize CBCT volumes using the same acquisition parameters as the clinical system. Offline CTsynthetic CBCT or synthetic CBCT -synthetic CBCT registration is then repeatedly performed while exhaustively searching for the optimal values of α and Γ. The search is then stopped when the configuration that provides the best similarity between the estimated displacements and the gold standard is found. An aspect which may prove to be slightly problematic for the EVolution algorithm is estimating motion in the proximity of sliding interfaces. The quadratic nature of the regularization term included in the functional from Eq. 1, constrains the estimated displacements to be spatially smooth. For soft tissues in particular, such an assumption is physically justified since they are incompressible and elastic and thus, shearing motion will in general not occur. On the other hand, at interfaces between organs and the abdominal or thoracic wall or between the organs themselves, sliding motion usually occurs, which locally violates the smoothness assumption. For example, during the experiment addressing CBCT-based respiratory motion compensation, this lead to a propagation of the displacements estimated for the kidneys and liver, ∼2 -3 voxels into the quasi-static abdominal and thoracic wall. However, this only had a local effect without notably affecting the overall performance of the algorithm. Nevertheless, as a source of misregistration this aspect may have to be monitored and can make the object of future studies. An important feature for registration algorithms in general is the required computational time. This is of particular interest when registration is employed for online therapy guidance, with the patient still on the interventional table.

Already from a numerical point-of-view, the EVolution algorithm facilitates an implementation with low computational demands since its optimization scheme is easily and massively parellelizable. This is due to the fact that the functional in Eq. 1 is an integral over a set of convex functions (one for each voxel of the image). It is a know result from convex optimization theory that the independent optimization of each of these functions leads to a global optimum of the functional. Therefore, the optimization of each function can and was delegated to a separate processing core, with all voxels being processed simultaneously. Moreover, while not illustrated in the presented work, all operations required by the numerical scheme of the EVolution algorithm can be implemented by using convolutions, which again render themselves easily parallelizable. As reported in section 3.6, for images of size 256 × 256 × 256 the algorithm requires ∼60 s per registration. This duration may scale with the image size, the extent of the displacements between the registered images and the available hardware. In order to correct for physiological motion and/or positioning errors, the current clinical standard for CBCT-guided radiotherapy implies the alignment of CT/CBCT images acquired at different stages of the treatment based on implanted fiducial markers. Besides causing discomfort to the patient during the implantation process, the marker-based alignment only accounts for rigid displacements, which may be suboptimal for deformable anatomies. By employing the EVolution algorithm instead, elastic deformation fields become available for both pathological tissues and organs-atrisk, providing the possibility of improved corrections due to the additional positional information. Moreover, this circumvents the implantation of fiducial markers, since the EVolution algorithm is purely image-based. We also anticipate that the dense elastic deformation fields could be used for adaptive replanning in the sense that the daily radiation dose could be delivered in accordance with the position of the pathology and the organs-at-risk at the time of treatment. Also, the daily dose can be projected with respect to the estimated displacements into the reference space of the planning CT image. This way, therapy progress could be monitored with respect to the initial plan and, if necessary, the remaining therapeutic dose can be adapted accordingly.

Conclusion

For an accurate and precise therapeutic dose delivery, image-guided radiotherapy treatments often rely on establishing a spatial coherence via image registration between CBCT and CT or CBCT images acquired at different stages of the therapy. Due to imaging-related radiation constraints, CBCT images are generally acquired using a low number of X-Ray projections and/or low beam intensities. This leads to acquisitions altered by streaking artifacts and a low SNR, which renders their registration using intensity-based methods challenging. The current study addresses both CT to CBCT and CBCT to CBCT registration by the means of the recently proposed EVolution algorithm. The method was validated in several complementary experiments where it was demonstrated that the artifacts associated to the sparse sampling of the CBCT sinogram only have a minimal impact, of sub-millimeter order, on the quality of the motion estimates. Compared to previous studies, the EVolution algorithm was proven to be capable of registering CBCT images of a considerably inferior quality. At the same time it did not require the integration of an additional intensity correction scheme. Moreover, the current study provides a robust quantification of the additional registration errors introduced by the artifacts present in the CBCT images. The EVolution algorithm requires the configuration of only two input parameters which, in practice, can be optimized prior to the therapy. Their values can then be maintained across different patients and anatomies, making the configuration of the method completely transparent to the clinician. In addition, the investigated method implies low computational demands, which is rarely the case for registration algorithms providing an elastic deformation between two images on a point-by-point basis. It can be therefore concluded that the EVolution algorithm is potentially an attractive asset for online motion compensation of CBCT-guided external beam radiotherapy treatments.

Figure 1 :

 1 Figure 1: Streaking artifacts in CBCT images: A transverse slice selected from a (a) CT (b) CBCT image acquired on a lung cancer patient. Compared to the CT image, the CBCT is visibly altered by streaking artifacts (indicated by the red arrows), which obstruct/degrade finer details of the anatomy and introduce intensity variations.

Figure 2 :

 2 Figure 2: (a) -(e) Effect of the number of projections on the CBCT images: Transverse slice from a CBCT image synthesized/reconstructed using (a) 360 (b) 120 (c) 72 (d) 50 (e) 36 projections. The original CT image used to synthesize the illustrated CBCT images, was acquired on a head-and-neck cancer patient. (f) -(j) Performance of the EVolution algorithm when registering a high resolution CT to CBCT images synthesized using a decreasing number of projections: The spatial distribution of the FEP between the estimated displacements and the gold standard, following the registration of a CT to a synthetic CBCT image reconstructed using (f) 360 (g) 120 (h) 72 (i) 50 (j) 36 projections. The alignment errors are illustrated as an overlay for the transverse slice in Fig. (a) -(e).

Figure 3 :

 3 Figure 3: Performance of the EVolution algorithm for registering a CT image from a head-and-neck cancer patient to a synthesized CBCT image: (a) -(c) A coronal, a sagittal and a transverse slice selected from the synthetic CBCT image. The red overlay indicates the location of the CTV. (d) -(f) The spatial distribution of the FEP before registration in the slices displayed in (a) -(c). (g) -(i) The spatial distribution of the FEP in the slices from (a) -(c) after registration. Note that for illustration purposes, two different color scales were used in the error maps from (d) -(f) and (g) -(i). Also, the patient's face was masked in the sagittal slice for privacy reasons.

Figure 4 :

 4 Figure 4: (a) -(b) A coronal, (e) -(f) a sagittal and (i) -(j) a transverse slice selected from a CBCT (first column) and a CT image (second column), acquired on a renal cell carcinoma patient. (c) -(d) A coronal, (g) -(h) a sagittal and (k) -(l) a transverse slice selected from a checkerboard overlap between the CBCT and the CT images before (third column) and after registration (fourth column).

Figure 5 :

 5 Figure 5: (a) -(b) A coronal, (e) -(f) a sagittal and (i) -(j) a transverse slice selected from a CBCT (first column) and a CT image (second column) acquired on a lung cancer patient. (c) -(d) A coronal, (g) -(h) a sagittal and (k) -(l) a transverse slice selected from a checkerboard overlap between the CBCT and the CT images, before (third column) and after registration (fourth column).

Figure 6 :

 6 Figure 6: Performance of the EVolution algorithm for CBCT-based respiratory motion estimation: (a) -(c) A coronal, a sagittal and a transverse slice selected from the CBCT image used as reference during the registration process. Point M indicates one of the implanted fiducial gold markers, while the blue contour delineates the contralateral kidney. (d) -(f) The spatial distribution of the temporally averaged FEP in the slices from (a) -(c), before aligning the synthetic 4D CBCT series to the reference image. (g) -(i) The spatial distribution of the temporally averaged FEP in the slices from (a) -(c), after registering the 4D series to the reference position.

Figure 7 :

 7 Figure 7: Performance of the EVolution algorithm for CBCT-based respiratory motion estimation: Trajectory of a fiducial gold marker and the centroid of the contralateral kidney over the respiratory cycle in the (a) -(b) head -foot, (c) -(d) anterior -posterior and (e) -(f) left -right directions. The red curves correspond to the trajectories estimated on the synthetic 4D CBCT series, while the blue curves plot their respective gold standards. The ordinate of the graphs indicate the displacement in millimeters, while the abscissa indicates the percentage phase within the respiratory cycle.

Figure 8 :

 8 Figure 8: Accuracy of the EVolution algorithm as a function of the input parameters α and Γ: (a) For CT -CT registration. (b) For CT -synthetic CBCT registration.

Figure 9 :

 9 Figure 9: Accuracy of the optical flow algorithm as a function of the regularization parameter β.

Table 1 :

 1 Table 2 also correspond to a larger DSC difference between the two methods in Table1. DSC, NCC and SSIM after the registration of pairs of CT images acquired on 5 head-and-neck cancer patients. The images were registered using both the optical flow (OF) and the EVolution (EVo) algorithm. Results are reported for several volumes-ofinterest (VOI). The missing values indicate that delineations were not available for the respective anatomies.

	Criterion	VOI	Method P#1 P#2 P#3 P#4 P#5
		Body	OF EVo	0.98 0.98 0.98 0.97 0.97 0.98 0.98 0.98 0.98 0.98
		Spinal Cord	OF EVo	0.77 0.8	--	0.84 0.7 0.78 0.83 0.71 0.8
		Right Parotid	OF EVo	0.85 0.85 0.85 0.71 0.73 0.86 0.84 0.86 0.81 0.63
	DSC	Left Parotid	OF EVo	0.84 0.83 0.88 0.74 0.77 0.85 0.82 0.9 0.82 0.76
		CTV	OF EVo	0.69 0.83 0.87 0.89 0.77 0.7 0.82 0.87 0.89 0.73
		Lymph Node #1	OF EVo	0.83 0.82 0.89 0.76 0.8 0.84 0.82 0.87 0.75 0.83
		Lymph Node #2	OF EVo	0.8 0.8	--	0.87 0.66 0.8 0.86 0.67 0.88
	NCC	Entire Volume	OF EVo	0.98 0.99 0.99 0.97 0.99 0.98 0.98 0.99 0.97 0.98
	SSIM	Entire Volume	OF EVo	0.93 0.96 0.97 0.92 0.93 0.92 0.95 0.96 0.91 0.92

Table 2 :

 2 FEP between the motion fields provided by the EVolution and the optical flow methods, respectively. Reporting is made under the format mean ± standard deviation, for the same CT data used to generate Table1.

Table 3 :

 3 Qualitative and quantitative evaluation of the algorithm performance when registering a CT to a CBCT image synthesized using 72 projections. Statistical distribution of the FEP is reported in columns 2 and 3 for the anatomical sites indicated in the first column, before and after registration. The fourth and fifth column of the table report the DSC for the indicated anatomies before and after registration.

	VOI	FEP before[mm] FEP after[mm] DSC DSC
		(mean±stdev)	(mean±stdev)	before after
	CTV	16.58 ± 2.2	0.31 ± 0.18	0.4	0.7
	Lymph node 1	15.85 ± 1.27	0.7 ± 0.37	0.27	0.84
	Lymph node 2	14.38 ± 1.68	0.65 ± 0.33	0.56	0.8
	Body	16.1 ± 3.45	0.66 ± 0.6	0.91	0.98
	Spinal cord	16.21 ± 1.74	0.15 ± 0.09	0	0.8
	Left parotid	16.38 ± 1.63	0.73 ± 0.42	0.38	0.84
	Right parotid	12.7 ± 1.12	0.65 ± 0.35	0.58	0.85

Acknowledgments

This work was supported by the Dutch Technology Foundation (STW) (project OnTrack #12813).