C Zachiu
email: c.zachiu@umcutrecht.nl

Baudouin Denis De Senneville

T Willigenburg

J Voort Van Zyp

J C J De Boer

B W Raaymakers

M Ries

J R N Voort Van Zyp

Anatomically-adaptive multi

Keywords: Multi-modal Deformable image registration, Anatomical adaptation, Imageguided radiotherapy

Introduction

One of the major challenges during external-beam radiotherapy (EBRT) is addressing the geometrical uncertainties introduced by the changes in shape and location of the tumor and the organs-at-risk (OARs) over the course of the treatment [START_REF] Roach | Diagnostic and therapeutic imaging for cancer: Therapeutic considerations and future directions[END_REF]. In case such uncertainties are not taken into consideration during the planning and delivery of EBRT, there is a high risk of under-dosage to the tumor, while at the same time over-irradiating adjacent healthy tissues [START_REF] Chavaudra | Définition des volumes en radiothérapie externe : rapports ICRU 50 et 62[END_REF][START_REF] Roach | Diagnostic and therapeutic imaging for cancer: Therapeutic considerations and future directions[END_REF][START_REF] Jaffray | Image-guided radiotherapy: from current concept to future perspectives[END_REF]). However, the recent integration of on-board imagers within the radiotherapy delivery systems has allowed visualizing the treated area and its surroundings during all phases of an EBRT work-flow: planning, delivery and response assessment [START_REF] Guckenberger | Image-guided radiotherapy based on kilovoltage cone-beam computed tomography -a review of technology and clinical outcome[END_REF], Raaymakers et al. 2017). This allows clinicians to identify the anatomical areas of interest on the acquired images and in turn, reduce the impact of geometric shifts and deformations on the overall treatment. Tracking the shape and location of the tumor and the OARs over the course of imageguided radiotherapy (IGRT) is typically done manually by experienced physicians [START_REF] Eisenhauer | New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1)[END_REF][START_REF] Mundt | Image-Guided Radiation Therapy: A Clinical Perpective[END_REF]. However, given the current tendencies towards imaging protocols which may include several imaging modalities and contrasts over the full course of the treatment, this manual process can become severely timeconsuming and error-prone. In addition, the time required for manual contouring may render a smooth clinical work-flow in an online setting unfeasible, in particular for onthe-fly correction strategies with the patient already on the treatment table. Therefore, an automatic tracking solution would be preferred instead. A feasible solution for automatic tracking of organ and pathological tissue boundaries over the course of the treatment is multi-modal deformable image registration (DIR) [START_REF] Hill | Medical image registration[END_REF][START_REF] Mani | Survey of Medical Image Registration[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. Such methods have the capability to estimate voxel-wise deformations across images acquired either with the same or a different modality and/or contrast. In effect, the patient anatomy can be practically tracked in an automatic manner within each acquired image, over the entire duration of the treatment. Moreover, the estimated deformations allow the up-stream mapping into the reference space of the therapy planning image of the dose delivered by each radiation fraction, allowing thus dose accumulation in a spatially consistent manner. This, in turn, gives way to potential adaptations of the therapeutic plan over the course of the therapy [START_REF] Kontaxis | Towards fast online intrafraction replanning for freebreathing stereotactic body radiation therapy with the MR-linac[END_REF]. A particularly attractive type of multi-modal DIR methods for IGRT are the socalled variational approaches [START_REF] Weickert | Variational optic flow computation: From continuous models to algorithms[END_REF]. Due to their high accuracy and precision, low number input parameters and rapid convergence, such algorithms were demonstrated to be especially beneficial in applications demanding on-the-fly corrections or even real-time performance (e.g. during therapy delivery), where short computational latencies are paramount [START_REF] Ries | Real-time 3d target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF][START_REF] Glitzner | On-line 3D Motion Estimation Using Low Resolution MRI[END_REF], Zachiu et al. 2017b[START_REF] Lafitte | Accelerating multimodal image registration using a supervoxel-based variational framework[END_REF], Zachiu et al. 2018). As a functioning principle, such methods estimate deformations between two or more images as the minimizer of a cost function comprising two terms: a data fidelity term and a regularization term. The data fidelity term quantifies the similarity between the images to be registered and decreases as the alignment of the images improves. Data fidelity terms for variational multi-modal DIR methods have been built around concepts such as mutual information [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF][START_REF] Maes | Medical image registration using mutual information[END_REF], modality-independent descriptors [START_REF] Heinrich | MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration[END_REF][START_REF] Reaungamornrat | MIND demons for MR-to-CT deformable image registration in image-guided spine surgery[END_REF]) and normalized gradient fields (Denis de [START_REF] Denis De Senneville | EVolution: An Edge -Based Variational Method for Non -Rigid Multi -Modal Image Registration[END_REF][START_REF] Spahr | Multimodal image registration for liver radioembolization planning and patient assessment[END_REF]. However, the minimization of the data fidelity term alone usually involves solving an under-determinate system of equations, which therefore leads to a non-unique solution or divergence of the algorithm. Therefore, the regularization term of the cost function adds further constraints on the system, rendering it over-determinate and implicitly the minimization problem becomes well-posed. From a practical perspective, the regularization term limits the degrees of freedom of the estimated deformations and, for anatomically meaningful deformations, should be chosen according to the deformation patterns and physical characteristics of the observed anatomy [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF]. Multiple previous studies have chosen spatial smoothness as a suitable constraint on the estimated deformations [START_REF] Viergever | A survey of medical image registration -under review[END_REF]. While in many cases this has proven to be a reasonable choice, particularly for aligning organ boundaries, it has been nevertheless demonstrated that within boundaries of elastic soft tissues this type of regularization may result in anatomically implausible deformations [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF][START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF]. More precisely, smoothness as the regularization constraint frequently leads to considerable compressions and expansions within the deformation vector field, which are physically implausible in such nearincompressible tissue structures. This, in turn, may have a direct impact on the precision and accuracy of the radiation dose mapping and accumulation process. In order to address this, previous studies have suggested replacing the smoothness constraint with a penalty on the Jacobian determinant of the deformations [START_REF] Rohlfing | Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint[END_REF][START_REF] Haber | Numerical methods for volume preserving image registration[END_REF][START_REF] Haber | Image registration with guaranteed displacement regularity[END_REF], Zachiu et al. 2018). While this has led to acceptable boundary alignment and more anatomically meaningful deformations, it does so only for near-incompressible anatomical regions. For areas which do undergo volumetric changes, incompressibility as a motion constraint frequently leads to misregistrations, in particular in purely gas/liquid-filled compartments of the human anatomy. Since variational algorithms evaluate their respective cost function over the entire field-of-view (FOV), the resulting misregistration of such problematic regions may affect the estimation accuracy in adjacent incompressible regions, or worse -should such a problematic region dominate the cost function -even prevent convergence entirely. In the current work we propose an anatomically-adaptive multi-modal variational DIR algorithm, which employs a local regularization depending on the anatomical properties of the observed anatomy. A test bench is also constructed, which demonstrates the benefits of such an approach compared to registration models which employ a global regularization. Moreover, despite the more complex nature of an anatomically adaptive approach, the numerical implementation has been optimized such that the performance remains entirely compatible with the requirements of online IGRT.

Method description

2.1. Investigated registration algorithms 2.1.1. DIR algorithms employing a global regularization In the scope of this work we have selected as a point-of-comparison to our anatomically-adaptive approach (which is described in detail in the subsequent paragraph), two previously-proposed algorithms which respectively also employ smoothness and incompressibility as a regularization constraint but, contrary to the new approach, assume validity of their respective regularization over the entire field-of-view (FOV). Both methods were previously validated for specific scenarios within the context of adaptive image-guided radiotherapy [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF][START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF]).

(i) The EVolution registration algorithm (abbreviated EVO in the scope of this manuscript), initially proposed in (Denis de [START_REF] Denis De Senneville | EVolution: An Edge -Based Variational Method for Non -Rigid Multi -Modal Image Registration[END_REF], estimates the deformation between two images I 1 and I 2 as the minimizer of the following functional:

E EV O ( u) = r∈Ω e -C( u( r)) + α ∇u 1 ( r) 2 2 + ∇u 2 ( r) 2 2 + ∇u 3 ( r) 2 2 (1)
where:

C(u( r)) = s∈Γ | ∇I 1 ( s) • ∇I 2 ( s + u( s))| s∈Γ ∇I 1 ( s) 2 ∇I 2 ( s + u( s)) 2 (2)
with u = (u 1 , u 2 , u 3 ) being the 3D displacement, Ω the image domain, r the voxel position, ∇ is the spatial gradient operator, • 2 is the Euclidean norm, α is a parameter linking the two terms of the functional and Γ is a symmetric cubic neighborhood around r. The data fidelity term (i.e. the first integrand) of the cost function in eq. 1 favors the alignment of strong gradients (e.g. edges) present in both I 1 and I 2 . The regularization term (i.e. the second integrand) addresses the ill-posedness of minimizing the data fidelity term alone, by assuming that the underlying anatomical deformations are smooth/differentiable. The regularization parameter α controls the amount of smoothness that the estimated deformations should showcase. Additional details on the numerical minimization of eq. 1 together with extensive validation of the algorithm can be found in (Denis de [START_REF] Denis De Senneville | EVolution: An Edge -Based Variational Method for Non -Rigid Multi -Modal Image Registration[END_REF]).

(ii) The incompressibility-regularized EVolution algorithm (referred to as EVI for the remainder of the manuscript), is a variation of EVO, initially designed as an improved solution for estimating anatomically meaningful deformations within elastic soft-tissues [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF]. Similar to EVO, it is a multi-modal variational image registration algorithm, with the deformations being estimated as the minimizer of:

E EV I ( u) = r∈Ω e -C( u( r)) + β (J( r + u( r)) -1) 2 (3) with J( r + u( r)) = 1 + ∂u 1 ( r) ∂x ∂u 1 ( r) ∂y ∂u 1 ( r) ∂z ∂u 2 ( r) ∂x 1 + ∂u 2 ( r) ∂y ∂u 2 ( r) ∂z ∂u 3 ( r) ∂x ∂u 3 ( r) ∂y 1 + ∂u 3 ( r) ∂z (4)
being the Jacobian determinant of the deformations. The rest of the terms maintain their meaning from eq. 2. The regularization term in eq. 3 constrains the estimated deformations to be incompressible, which has been demonstrated to be a better choice within near-incompressible anatomical structures compared to EVO [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF][START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF]).

Proposed anatomically adaptive DIR solution

In the scope of this paper we propose replacing the global regularization employed by EVO and EVI, with a localized regularization which depends on the specific anatomical and physical properties of the observed anatomical structures. In effect, this new registration model which we will call adaptive-EVolution (AEVO) provides the deformations between two images as the minimizer of the following functional:

E AEV O ( u) = r∈Ω e -C( u( r)) + M S ( r) • γ ∇u 1 ( r) 2 2 + ∇u 2 ( r) 2 2 + ∇u 3 ( r) 2 2 + M I ( r) • δ (J( r + u( r)) -1) 2 (5)
where M S and M I are binary 3D matrices, defining the regions where the observed anatomy undergoes smooth or incompressible deformations, respectively. Similar to eq. 2 and 3, the regularization parameters γ and δ control the amount of smoothness and incompressibility within their respective image regions. The two binary masks M S and M I were chosen here to be complementary (i.e. M I = 1 -M S ). Alternatively, both regularizations could be employed simultaneously (i.e. M S = M I = 1) and allow parameters γ and δ to be spatially variant instead. The amount of smoothness and incompressibility of the estimated deformations could then be weighted in accordance with the biomechanical properties of the observed tissues.

The minimization of the cost function in eq. 5 was performed by running an iterative fixed-point scheme on the resulting Euler-Largrange equations, similar to the approach described in (Denis de [START_REF] Denis De Senneville | EVolution: An Edge -Based Variational Method for Non -Rigid Multi -Modal Image Registration[END_REF]) and [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF] for EVO and EVI. The registration process was performed in a coarse-to-fine manner, by iterating it from a 16-fold downsampled version of the images up to their original resolution, with an upsampling factor of two. Together with the images, the binary masks M S and M I were also downsampled accordingly, using and averaging downsampling kernel. At each resolution level, we have also employed iterative refinement of the deformations, which implied restarting the registration process several times at the same resolution, using as an initial value the deformations from the previous refinement iteration. The registration process at each resolution level was stopped when the relative average difference between the motion fields estimated by two successive iterations was lower than 5% (chosen empirically). The stopping criterion was evaluated individually for both the compressible and incompressible areas. Note that the EVO and EVI cost functions were minimized using the approach described in their original manuscripts (Denis de [START_REF] Denis De Senneville | EVolution: An Edge -Based Variational Method for Non -Rigid Multi -Modal Image Registration[END_REF], Zachiu et al. 2018).

Experimental setup

The AEVO, EVO and EVI registration methods were evaluated within the context of two clinically relevant applications of IGRT: CT/conebeam CT(CBCT)-guided radiotherapy for lung cancer and MR-guided radiotherapy of prostate cancer. The two anatomical sites were also chosen due to the considerably different deformation patterns showcased by the respective structures present within the field-of-view, i.e. compressible areas (lungs, bladder and rectum) situated in the immediate vicinity of near-incompressible structures (such as the pathological tissue itself). In the following we will describe the medical image data used for this purpose, together with the employed evaluation methodology for the algorithms.

Selection of the medical image datasets

For CT/CBCT-guided radiotherapy for lung cancer, we have selected a 4D-CT image sequence from five patients, with each sequence sampling the 3D anatomical changes over the course of their respiratory cycle (from 0 to 90 %, with an increment of 10%). The data was downloaded from "The Cancer Imaging Archive", originally collected in the scope of the following studies [START_REF] Roman | Interfractional positional variability of fiducial markers and primary tumors in locally advanced non-small-cell lung cancer during audiovisual biofeedback radiotherapy[END_REF][START_REF] Clark | The cancer imaging archive (TCIA): Maintaining and operating a public information repository[END_REF][START_REF] Balik | Evaluation of 4-dimensional computed tomography to 4-dimensional cone-beam computed tomography deformable image registration for lung cancer adaptive radiation therapy[END_REF][START_REF] Hugo | Data from 4d lung imaging of nsclc patients[END_REF][START_REF] Hugo | A longitudinal four-dimensional computed tomography and cone beam computed tomography dataset for image-guided radiation therapy research in lung cancer[END_REF]. For each dynamic image volume of the selected 4D-CT sequences, the curators of the database have also provided contours for the lungs, the tumor and several other structures. Out of each 4D-CT sequence we have only selected image pairs which showcase the largest anatomical differences between them. In effect, the fullinhalation image (the 0% phase) was established as reference image for the registration process, whereas the 40 -70 % images were considered to be due for registration. The resulting dataset was then used to evaluate the performance and behavior of EVO, EVI and AEVO for CT to CT registration. Note that for the AEVO method, the lungs, the heart and the area outside the body were considered to be compressible, while the rest of the structures present in the FOV were deemed near-incompressible.

In order to evaluate the performance of the algorithms for CT to CBCT registration (potentially beneficial for applications described in detail in (Zachiu et al. 2017b)), a set of CBCT images were synthesized from the 40 -70% CT images mentioned above. This was achieved by a three-fold down-sampling of their sinogram and subsequent reconstruction via the TIGRE toolbox [START_REF] Biguri | TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction[END_REF]. This resulted in a set of synthetic CBCT images with a quality similar to that provided by a regular CBCT imaging device (see (Zachiu et al. 2017b) for further discussion on this topic). Thus, registering the synthetic 40 -70% CBCT images to the original 0% CT via the three algorithms, allowed their evaluation for CT-CBCT registration. In terms of acquisition parameters, the original CT images were acquired on a Varian Trilogy scanner, in-plane image size 512 × 512 and a variable number of slices, voxel size 0.97 × 0.97 × 3.0mm 3 . The evaluation of the algorithms for MR-guided radiotherapy for prostate patients was performed on image data acquired as part of a cohort study on the Unity MR-Linac installed at the UMC Utrecht, Utrecht, The Netherlands. In-line with the potential requirements of adaptive IGRT on the MR-Linac, we have evaluated the algorithms for: 1) Mono-modal MR image registration; 2) Multi-modal/cross-contrast MR image registration and 3) CT to MR image registration. For each of these instances, we have selected image pairs from five patients, with contours for the prostate, bladder, rectum, femurs and the pelvis being available for each individual image. The contours were validated for accuracy by two experienced radiation oncologists. For registration using the proposed AEVO method, the prostate and the bony structures (i.e. the pelvis and the femurs) were considered to be near-incompressible, while the remainder of the areas (including the region outside the body) were allowed to smoothly deform. The acquisition parameters for the reference and the moving images in the three datasets are reported in had a different FOV and/or voxel size, the moving image was initially re-sampled on the grid of the reference. This was subsequently followed by resizing both images to 256 × 256 × 256, for computational purposes. For the CT -MR dataset, we have also performed an approximate rigid manual re-alignment of the images prior to applying deformable registration.

2.2.2.

Algorithm evaluation criteria The three algorithms were comparatively evaluated against three criteria: contour propagation performance, anatomical plausibility of the estimated deformations and accuracy and precision of the estimated deformations.

(i) In order to determine the performance of the algorithms for contour propagation, we employed the dice similarity coefficient (DSC) before and after registration [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF]:

DSC(A, B) = 2 * |A ∩ B| |A| + |B| (6)
where A and B are two contours whose overlap is determined and | • | corresponds to the number of voxels within a contour.

(ii) For anatomical plausibility, the voxel-wise Jacobian determinant of the estimated deformations was evaluated.

It is known from continuum mechanics that deformations of near-incompressible materials have a Jacobian determinant that is close to one. Similarly, within anatomical areas which are known to undergo volumetric changes (such as the lung, bladder and rectum), a value which is close to one may be indicative of misregistrations. For the mathematical expression of the Jacobian determinant of a 3D deformation field we refer the reader to eq. 4.

(iii) For the comparative evaluation of the accuracy and precision of the algorithms, we have applied a set of known deformations on one of the CT lung cancer images and one of the MR prostate cancer images. In this manner, two pairs of synthetically deformed images were generated and subsequently registered to one-another via the investigated algorithms. The voxel-wise root square error between the known and the estimated deformations was then calculated and used as an evaluation criterion:

RSE( r) = u est ( r) -u known ( r) 2 (7)
where u est and u known are the estimated and the known deformations, respectively,

• 2 is the Euclidean norm and r is the voxel position. The known deformations were generated via a finite element modeling (FEM) of the organ displacements in each of the two images. In summary, the available organ contours were initially imported into the FEM software PreView v2.1.4 [START_REF] Maas | FEBio: Finite elements for biomechanics[END_REF], which allowed establishing the physical properties of the tissues and the motion actuators. The latter implied a lung inflation for the CT data and an increase of the bladder volume and rectum for the MRI. This was followed by the effective finite element simulation of the resulting deformations via the FEBio v2.9.1 software [START_REF] Maas | FEBio: Finite elements for biomechanics[END_REF]. Finally, the simulated deformations were inspected using PostView v2.4.4 [START_REF] Maas | FEBio: Finite elements for biomechanics[END_REF]) in order to confirm their similarity to the typical deformation patterns occurring in these anatomical regions. The resulting vector fields were then used to deform the original CT and MR image, respectively. Note that the employed elastic modulus (E) and the Poisson ratio (ν) for the tissues of interest, required as part of the simulation, were adopted from previous studies (Qiao et

Algorithm implementation and calibration

All algorithms were implemented using the Compute Unified Device Architecture (CUDA) and executed on a Nvidia TITAN V graphics card. The regularization parameters α, β, γ and δ (see eq. 1, 3 and 5) were configured by an exhaustive search, maximizing the post-registration contour alignment of the clinical target volume (CTV). The calibration procedure was only performed on one of the image pairs in each of the datasets (CT-CT, CT-CBCT, MR Monomodal, MR Multimodal and CT-MR), with the parameter values being maintained for the remainder of the image pairs in each dataset. Such an optimization procedure is justified by the fact that in clinical practice, registration results are validated by a visual inspection of the post-registration organ boundary alignment. In effect, we have used the maximization of the post-registration DSC as an objective surrogate for such a criterion. The employed values for the regularization parameters are reported in table 2. The size of the neighborhood Γ required for the calculation of the data fidelity term was set to 11 × 11 × 11 for all algorithms, as suggested in (Denis de [START_REF] Denis De Senneville | EVolution: An Edge -Based Variational Method for Non -Rigid Multi -Modal Image Registration[END_REF]) and [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF]. Table 2: Regularization parameters for the three registration algorithms, for each of the image modalities.

Dataset

Results

3.1. Performance of the investigated algorithms for tracking organ boundaries 3.1.1. CT/CBCT-based tracking in lung cancer patients Figure 1 illustrates the statistical distribution of the Dice Similarity Coefficient (DSC) following the registration of the CT-CT and CT-CBCT datasets (see section 2.2.1). The pre-and post-registration DSC values were pooled together from the five patients, for each of the anatomical structures-of-interest, and displayed under the shape of a boxplot. The box limits of each boxplot are represented by the 25th and 75th percentiles, while the whiskers are approximately the 1st and the 99th percentiles, respectively. While all three algorithms lead to improvements of the post-registration DSC, there is a noticeable tendency of the EVI algorithm to under-perform for the lungs, compared to both EVO and AEVO.

On the other hand, there are no considerable differences between the algorithms for the tumor itself (assumed near-incompressible). The two images showcase the statistical distribution of the DSC for the lungs and the tumor before and after registration using the three algorithms.

MR-based tracking in prostate cancer patients

2.2.1).

The values in the table represent the DSC averaged over the five image pairs from each of the MR Monomodal, MR Multimodal and CT-MR datasets. The DSC values are reported individually for the bladder, prostate and rectum. Similar to the CT/CBCT tracking scenario, all three algorithms have led to a notable improvement of the DSC in all cases, compared to pre-registration. On the other hand, while for the prostate itself all three algorithms perform similarly, for the bladder and rectum, the EVI algorithm has a tendency to misregister the images in these areas, with differences of the average DSC of up to 7% compared to EVO and AEVO. The latter two algorithms, however, showcase no notable differences in terms of the DSC.

ROI

Anatomical plausibility of the estimated deformations

Figure 2 illustrates the voxel-wise Jacobian determinant of the deformations estimated by the evaluated algorithms on one of the image pairs from each of the five datasets. It can be observed that in all cases, the EVO algorithm displays moderate to high deviations of the Jacobian determinant from one, regardless whether the underlying anatomy is near-incompressible or not. On the other hand, as the EVI algorithm penalizes deviations of the Jacobian determinant from one, such deviations are considerably reduced over the entire field-of-view. This also includes anatomical structures which are expected to change volume such as the lung and the rectum.

The AEVO method, on the other hand, maintains a Jacobian determinant close to one for structures which were considered to be near-incompressible in the scope of this study (enclosed by the red dashed lines), while for the remainder of the structures moderate to high deviations can be observed. Another interesting observation is that for the illustrated CT-CBCT case, the AEVO algorithm showcases a more uniform spatial distribution of the Jacobian determinant within the lung, compared to the EVO method.

To support the observations stemming from figure 2, figure 3 shows the statistical distribution of the Jacobian determinant provided within incompressible structures by the three registration algorithms. The illustrated boxplots were generated by pooling together the Jacobian determinant of the deformations within incompressible structures from all registered image pairs in each of the five datasets. It can be observed that, to different extents, the EVO algorithm leads in all cases to considerably larger variations of the Jacobian determinant compared to EVI and AEVO. On the other hand, within incompressible structures, AEVO demonstrates similar ranges of the Jacobian determinant to EVI.

Algorithm performance with respect to known deformations

Figure 4 showcases a selected slice from the anatomical images used in the FEM experiment (see section 2.2.2). More precisely, the images illustrate a coronal (figures 4(a)-4(c)) and a sagittal (figures 4(d)-4(f)) slice selected from the reference image, the moving image and their color-coded fusion, associated to the lung CT and prostate MR FEM simulation. Within the color-coded fusion, the magenta channel corresponds to the moving image, while the green channel is the reference. As intended, it can be observed that the anatomical differences between the images stem from an inflation of the lungs for the CT data and a volumetric increase of the bladder and rectum for the prostate cancer data. The spatial distribution of the Jacobian determinant of the deformations estimated by the three algorithms on the FEM-generated data is illustrated in figure 5. The results are consistent with the observations made for the clinical data: the EVO and AEVO algorithms showcase high deviations from one within compressible structures such as the lungs, bladder and rectum, whereas the values for EVI remain in the proximity of one. On the other hand, within incompressible anatomies, the EVO method showcases moderate to high deviations from one, while AEVO and EVI provide similar values. This is again re-confirmed by the statistical analysis of the Jacobian determinant within incompressible structures, showcased within figure 6.

Figure 7 displays the statistical distribution of the RSE calculated between the estimated and the simulated deformations, for all three algorithms. The errors were aggregated from both the compressible and incompressible structures in order to generate the illustrated boxplots. It can be noted that for both image pairs, the AEVO algorithm has a tendency to provide lower estimation errors compared to both EVO and EVI, with a maximum RSE of 1.5 -2 mm. On the other hand, the EVI method leads to rather large errors of up to ∼7 mm for the CT data and ∼9 mm for the MR images. The EVO algorithm is situated in-between AEVO and EVI, with errors of up to ∼4.5 mm for the lung images and ∼2.5 mm for the prostate data, with the RSE for the latter being overall only slightly higher compared to AEVO.

Computational time of the proposed registration algorithm

Table 4 reports the average computational time required by the proposed AEVO registration algorithm, for each of the five registered image modalities. Recall that all images were re-sampled on a 256 × 256 × 256 grid prior to registration (see section 2.2.1). It can be observed that in all instances, the average convergence time remains under one minute. Noteworthy is also the fact that for CT-to-MR image registration, the algorithm converged considerably faster to a solution, compared to the rest of the imaging modalities.

Discussion

Several applications within image-guided radiotherapy such as daily positioning, contour propagation, dose accumulation and adaptive re-planning may potentially benefit from the inclusion of DIR algorithms in their workflows, with the aim of improving the geometric accuracy of treatment planning, delivery and response assessment [START_REF] Brock | Use of image registration and fusion algorithms and techniques in radiotherapy: Report of the AAPM radiation therapy committee task group no. 132[END_REF]. Since current IGRT may include multiple imaging modalities and contrasts over the course of the treatment, multi-modal DIR methods are of particular interest in this context. Additionally, in case such methods are employed for applications requiring on-fly-processing, for example if the patient is already positioned on/in the therapy system, computational times and robust convergence become an important factor. In this sense, variational multi-modal DIR algorithms are a particularly attractive solution, due to their previously demonstrated precision, accuracy and fast convergence times Here, we propose a variational multi-modal registration algorithm, which adapts its registration model in accordance with the observed underlying anatomy (AEVO). In order to determine the benefits such an approach may provide, we have comparatively evaluated it against two existing solutions: one with a regularization term stemming from the field of digital image processing (EVO) (Denis de Senneville et al. 2016) and its variation which has been adapted for tracking elastic soft tissues (EVI) [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF]). The rationale is hereby that the majority of the existing cross-modality variational registration algorithms such as EVO or EVI employ a spatially invariant regularization across the entire field of view. The underlying assumption is hereby, that (a) either the overall tissue structure is sufficiently homogenous with respect to their biomechanical properties to justify this generalized choice of the regularization, or (b) that the data fidelity term of the variational is strong enough "to pull" the algorithm even in those regions to the anatomically correct result, where assumption (a) is (partially) violated. As our results and the results of numerous papers have shown, this assumption is generally well fulfilled for the evaluated EVO and EVI algorithms and can furthermore be addressed by choosing the appropriate regularization parameter. However, in the direct vicinity of biomechanically heterogeneous regions, and in particular for imaging modalities (or anatomical regions) associated with poorer soft tissue contrast, these assumptions become problematic due to the limited contribution of the data fidelity term to the variational in these regions. As it has been shown in previous studies for finite element-based registration algorithms [START_REF] Bharatha | Evaluation of three-dimensional finite element-based deformable registration of pre-and intraoperative prostate imaging[END_REF][START_REF] Brock | Accuracy of finite element model-based multi-organ deformable image registration[END_REF][START_REF] Zhong | A finite element method to correct deformable image registration errors in low-contrast regions[END_REF][START_REF] Velec | Validation of biomechanical deformable image registration in the abdomen, thorax, and pelvis in a commercial radiotherapy treatment planning system[END_REF], taking local biomechanical properties as prior knowledge into account can significantly improve the result of the registration process. Unfortunately, for most biomedical image modalities such as MRI, CT or CBCT it is generally not possible to derive precise local knowledge of all biomechanical parameters on an inter-individual basis from the images themselves. Nevertheless, almost all these imaging modalities allow a coarse classification of imaged tissue into the categories air, liquid, soft-tissue and bone. This, in-turn, unlocks the possibility to employ variational algorithms with a locally variant regularization, which matches locally best for a given type of tissue.

An additional point to be made is that all three algorithms employ the same data fidelity term, allowing thus to study the impact the different regularization terms have on the estimated deformations.

Algorithm performance for contour propagation

The contour propagation capabilities of all three algorithms was evaluated in the context of registering multiple imaging modalities which may be involved in an IGRT workflow: CT -CT, CT -CBCT, same-contrast MR, multi-contrast MR and CT -MR. The images were acquired as part of a protocol for lung and prostate cancer patients, respectively. Both of these anatomical sites contain neighboring regions with considerably different deformation patterns. For example, while the lungs can undergo significant volumetric changes over the respiratory cycle, the surrounding tissues are only displaced/transported from one part of the FOV to another, with limited changes in volume. A similar situation arises in the prostate, which may be displaced under the effect of volumetric changes within the bladder and/or rectum.

In order to determine the contour propagation performance of the three algorithms, we have evaluated the pre-and post-registration dice similarity coefficient (DSC) for several structures of interest in each of the five datasets. As shown in figure 1 andtable 3, all algorithms have led to a notable improvement of the DSC compared to the preregistration case. However, for anatomies which can undergo volumetric changes (such as the lungs, bladder and rectum), the EVI algorithm leads to systematically larger errors in terms of contour alignment, compared to EVO and AEVO. Such an outcome is, nevertheless, in line with the design of EVI, which globally penalizes deviations of the Jacobian determinant from one and implicitly any volumetric changes. Concerning AEVO, while there is a tendency to provide slightly lower values for the DSC compared to EVO, this situation only occurs within isolated cases, with limited differences of 1-2 %.

On the other hand, for near-incompressible structures such as the pathological area, all three algorithms provided similar DSC values in both the lung and the prostate cancer cases. It is worth noting that for EVO and EVI this is in good correspondence with our previous findings in [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF][START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF].

Anatomical plausibility of the estimated deformations

As previously discussed, the Jacobian determinant of near-incompressible anatomies should be close to one. Thus, large deviations from this value in such areas are physically implausible and indicative of misregistrations [START_REF] Schreibmann | A measure to evaluate deformable registration fields in clinical settings[END_REF], Zachiu et al. 2018). In turn, this may lead to mapping errors of quantitative information such as Houndsfield units, radiation dose and/or diffusion/perfusion values, as discussed and demonstrated by previous studies [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF][START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF]. In effect, we have evaluated here the spatial and statistical distribution of the Jacobian determinant for all three of the algorithms. What is already apparent from a visual analysis of figures 2 and 5 is that both EVO and EVI provide a spatial distribution of the Jacobian determinant with consistent characteristics over the entire FOV. While EVO showcases moderate to high deviations from one within both compressible and near-incompressible structures, EVI overall penalizes such deviations, including within structures which typically undergo volumetric changes. Particularly for the FEM experiments (see figure 5), where large volumetric changes were simulated for the lungs and the bladder, a Jacobian determinant close to one within these structures is a clear indication of misregistration. AEVO, on the other hand, limits such deviations solely within structures which were considered to be near-incompressible, with the areas outside these structures showcasing volumetric changes to approximately the same extent as EVO. A noteworthy observation is that, within "compressible" regions, the AEVO algorithm leads to smoother spatial variations of the Jacobian determinant compared to EVO. This effect is particularly observable within the lungs and is in good correspondence with the physical reality, since during respiration the pressure within the lungs is near-uniform. We hypothesize that this stems from the multi-resolution scheme used to minimize the AEVO variational (see section 2.1.2). At the lower resolution levels, both the smoothness and incompressible regularizations may end-up operating on the same voxels, particularly for the ones situated at the boundary between the two regions. This can lead to an incompressible deformation pattern, partially being present within smooth regions of the motion field used as an initializer for the higher resolutions. In turn, this may lead to a more uniform spatial distribution of the Jacobian determinant within compressible areas.

The above observations on the Jacobian determinant values for the three algorithms were confirmed by its statistical analysis within near-incompressible regions. It can be seen in figures 3 and 6 that while AEVO and EVI provide values within a similar range, EVO leads to a systematic ∼20 -40% higher range of the Jacobian determinant. In effect, the deformations provided by AEVO and EVI are more anatomically plausible in near-incompressible anatomical structures, compared to EVO.

Overall performance of the algorithms

Depending on the specific particularities and requirements of the application at hand, a registration algorithm may or may not be suitable for the task. In case contour propagation within elastic soft tissues is of interest, the results in figure 1 and table 3 indicate that either of the three evaluate algorithms would provide satisfactory results.

For anatomical structures which undergo volumetric changes however, the EVI method is by design unsuitable to handle such deformation patterns.

On the other hand, in case the application requires the mapping of quantitative information within elastic soft tissues, the EVO algorithm may lead to mapping errors due to its tendency to provide moderate to high deviations of the Jacobian determinant from one [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF][START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF], while AEVO and EVI considerably mitigate this effect (see figures 3 and 6). Nevertheless, within compressible anatomical structures, EVI is unsuitable for this task as well, due to its limited ability to register volumetric changes, in particular within structures displaying poor soft-tissue contrast.

In effect, the proposed anatomically-adaptive approach appears to be the optimal choice among the three: it is able to accurately track the boundaries of both compressible and incompressible anatomical areas, while at the same time maintaining the anatomical plausibility of the estimated deformations within elastic soft tissues. This assertion is also supported by the two FEM experiments, which have demonstrated that in both cases the AEVO algorithm outperforms in both accuracy and precision (to different extents) the EVO and EVI methods (as shown in figure 7).

Compared to EVO and EVI, however, the AEVO algorithm requires a map of compressible/incompressible regions as an input. This basically provides the registration method with prior information regarding the deformation characteristics of different anatomical regions present within the FOV. While this may seem as a limitation of the approach, contours of the tumor and the organs-at-risk (allowing the definition of such areas) are readily available from the therapy planning step of the IGRT. Additionally, the recent developments within the domain of deep learning-based segmentation, present the potential to automate the generation of the required maps.

Another distinguishable difference between AEVO and the two other methods is the requirement to calibrate two regularization parameters instead of one (see eq. 1, 3 and 5). While this leads to longer calibration times, we have observed that from a practical point-of-view, this only needs to be performed once per registration modality. In effect, similar to EVO and EVI, the regularization parameters can be prospectively optimized and, as long as the MR/CT acquisition parameters do not change significantly, the same values can be maintained over practically any number of cases. Nevertheless, for future studies we intend to investigate methodologies which allow the automatic calibration of the proposed AEVO method.

In terms of computational requirements, the average convergence time for AEVO remained under 60 s for all registered image modalities. An interesting observation is that CT-to-MR registration was considerably faster compared to the rest of the imaging modalities (see table 4). We hypothesize that this is due to the strong dissimilarity in the appearance of anatomical structures between the CT and the MR images. In effect, it is only the most prominent features such as organ boundaries which contribute to the data fidelity term of the AEVO cost function in eq. 5, requiring overall less iterations to optimize. Nevertheless, the convergence time remains compatible with clinical workflows in an online setting. Naturally, the computational requirements may scale with the size of the images, the employed hardware and magnitude of the deformations.

Conclusion

The present study proposes a novel multi-modal variational DIR algorithm, which adapts its registration model in accordance with the local properties of the observed anatomy. The method was comparatively evaluated against two existing DIR methods EVO and EVI, which employ a smoothness and an incompressible regularization over the entire field-of-view, respectively. In terms of contour propagation, both the proposed AEVO and the existing EVO have shown a similar performance, whereas EVI was demonstrated to be sub-optimal for areas which undergo volumetric changes.

On the other hand, while the EVO algorithm can lead to anatomically implausible deformations within elastic biological soft-tissues, the proposed AEVO method addresses this shortcoming, leading to a local Jacobian determinant similar to that of EVI.

As demonstrated in previous studies, this leads to a more precise and accurate mapping of quantitative information such as the delivered radiation dose. The FEM experiments have also demonstrated that overall, the AEVO algorithm has the potential to estimate deformations with a higher accuracy and precision, compared to the existing methods. Moreover, despite the more complex registration model and numerical implementation, the proposed anatomically adaptive approach leads to computational times which remain compatible with clinical scenarios in an online setting. Therefore, we can conclude that, compared to the methods employing a global regularization, the proposed AEVO algorithm showcases better potential benefits for future adaptive IGRT workflows.

  al. 2005, Chai et al. 2009, Al-Mayah et al. 2009, Lee et al. 2012, Zachiu et al. 2020).

Figure 1 :

 1 Figure 1: Boundary tracking capabilities of the investigated methods for (a) CT-based and (b) CBCT-based tracking. The two images showcase the statistical distribution of the DSC for the lungs and the tumor before and after registration using the three algorithms.

Figure 2 :

 2 Figure 2: Anatomical plausibility of the deformations estimated by the three algorithms on the five datasets -example. The first column showcases a slice from the reference image of one of the registered image pairs from each of the five datasets. Columns two, three and four display the spatial distribution of the Jacobian determinant of the deformations estimated by EVO, EVI and AEVO, respectively. The red dashed lines indicate the anatomical areas which were considered to be incompressible in the scope of this work. The blue dashed lines within the first two rows enclose the lungs, which were considered compressible.

Figure 3 :

 3 Figure 3: Statistical distribution of the Jacobian determinant provided by the three algorithms within incompressible areas. Illustration is performed individually for each of the five datasets.

Figure 4 :

 4 Figure 4: Anatomical images employed in the FEM experiment. The figures illustrate a selected slice from the reference, moving and color-fused image for (a)-(c): the lung CT and (d)-(f): the prostate MR FEM simulation.

Figure 5 :Figure 6 :

 56 Figure5: The spatial distribution of the Jacobian determinant for each of the three analyzed algorithms when applied on the FEM-generated image pairs. The illustration is made for the same slices as in figure4. The red dashed lines enclose the anatomical structures which were considered to be incompressible, i.e. the tumor, rib cage, liver and spleen for the CT data, while for the MR data it is the prostate and the pelvis. The blue dashed line indicates the compressible structures, namely the lungs and the heart for the CT image and with the bladder and the rectum for the MR data.

(Figure 7 :

 7 Figure 7: Statistical distribution of the RSE between the deformations estimated by the three algorithms and the FEM-simulated deformations.

Table 1 :

 1 Acquisition parameters for the images used in the scope of evaluating EVO, EVI and AEVO for MR-guided radiotherapy of prostate cancer patients. The employed abbreviations are as follows: TR -repetition time, TE -echo time and FA -flip angle.

	Dataset	Image	TR [ms] [ms] TE	FA •	Image Dimension	Voxel Size [mm 3 ]
	MR-MR Reference 1535	277	90	480 × 480× 300 0.83 × 0.83 × 2
	Mono	Moving			Same as the reference image
	MR-MR Reference 1535 277.81	90	480 × 480× 300 0.83 × 0.83 × 2
	Multi	Moving	4.68	2.34	50	448 × 448× 63 0.93 × 0.93 × 2
	CT-MR	Reference Moving	-1535	-277	-90	512 × 512× 125 1.03 × 1.03 × 3 480 × 480× 300 0.83 × 0.83 × 2

table

1

. Note that all MR images were acquired on a 1.5 T Philips Achieva scanner, Philips Healthcare, Best, The Netherlands, while the CT prostate cancer images were acquired on a Philips Brilliance Big Bore scanner, Philips Healthcare, Best, The Netherlands. For instances in which the reference and the moving image

  Table 3 reports the DSC before and after registration of the datasets acquired on the prostate cancer patients (see section

		1
		0.95
	Dice Similarity Coefficient	0.8 0.85 0.9
		0.75
		0.7
		Right Lung	Left Lung	Tumor

None EVO EVI AEVO None EVO EVI AEVO None EVO EVI AEVO

Table 3 :

 3 Boundary tracking capabilities of the evaluated algorithms for MR-based guidance in the context of EBRT for prostate cancer. The values represent the average DSC over the five image pairs in each of the MR Monomodal, MR Multimodal and CT-MR datasets. The reporting is made individually for the bladder, prostate and rectum.

		Algorithm MR-MR Mono MR-MR Multi MR -CT
		None	0.79	0.79	0.73
	Bladder	EVO EVI	0.93 0.86	0.91 0.85	0.86 0.81
		AEVO	0.93	0.90	0.85
		None	0.83	0.78	0.78
	Prostate (CTV)	EVO EVI	0.92 0.92	0.83 0.84	0.82 0.83
		AEVO	0.93	0.84	0.82
		None	0.75	0.79	0.73
	Rectum	EVO EVI	0.88 0.82	0.85 0.82	0.80 0.80
		AEVO	0.88	0.85	0.80

Acknowledgments

This work was supported by ITEA 3, project no. 16016 (STARLIT).