Igor Peterlík

Hadrien Courtecuisse

Robert Rohling

Purang Abolmaesumi

Christopher Nguan

Stéphane Cotin

Septimiu Salcudean

Fast Elastic Registration of Soft Tissues under Large Deformations

A fast and accurate fusion of intra-operative images with a pre-operative data is a key component of computer-aided interventions which aim at improving the outcomes of the intervention while reducing the patient's discomfort. In this paper, we focus on the problematic of the intra-operative navigation during abdominal surgery, which requires an accurate registration of tissues undergoing large deformations. Such a scenario occurs in the case of partial hepatectomy: to facilitate the access to the pathology, e.g. a tumor located in the posterior part of the right lobe, the surgery is performed on a patient in lateral position. Due to the change in patient's position, the resection plan based on the pre-operative CT scan acquired in the supine position must be updated to account for the deformations. We suppose that an imaging modality, such as the cone-beam CT, provides the information about the intra-operative shape of an organ, however, due to the reduced radiation dose and contrast, the actual locations of the internal structures necessary to update the planning are not available.

To this end, we propose a method allowing for fast registration of the pre-operative data represented by a detailed 3D model of the liver and its internal structure and the actual configuration given by the organ surface extracted from the intra-operative image. The algorithm behind the method combines the iterative closest point technique with a biomechanical model based on a co-rotational formulation of linear elasticity which accounts for large deformations of the tissue. The performance, robustness and accuracy of the method is quantitatively assessed on a control semi-synthetic dataset with known ground truth and a real dataset composed of nine pairs of abdominal CT scans acquired in supine and flank positions. It is shown that the proposed surface-matching method is capable of reducing the target registration error evaluated of the internal structures of the organ from more than 40 mm to less then 10 mm. Moreover, the control data is used to demonstrate the compatibility of the method with intra-operative clinical scenario, while the real datasets are utilized to study the impact of parametrization on the accuracy of the method. The method is also compared to a state-of-the art intensity-based registration technique in terms of accuracy and performance.

Introduction

A vast amount of research has been dedicated to deformable registration due to its potential clinical impact. Many methods have been proposed to address this complex task, and the solutions usually depend on the actual scenario and the type of data to be registered [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. Image registration has been widely addressed in the context of radiation therapy, intervention planning, intra-operative navigation and others.

The computer-based methods related to patient diagnosis and pre-and post-operative treatment represent an important set of applications requiring accurate and reliable registration of medical images [START_REF] Hill | Medical image registration[END_REF][START_REF] Khalifa | State-of-the-art medical image registration methodologies: a survey[END_REF][START_REF] Crum | Non-rigid image registration: theory and practice[END_REF]. In the context of computer-aided diagnosis, the medical image registration is applied to two or more images acquired at different times in order to quantify the evolution of tissue pathology or response to treatment [START_REF] Charnoz | Liver registration for the follow-up of hepatic tumors[END_REF][START_REF] Kaus | Assessment of a model-based deformable image registration approach for radiation therapy planning[END_REF][START_REF] Eisenhauer | New response evaluation criteria in solid tumours: revised recist guideline (version 1.1)[END_REF]. In computed-aided plan-ning and navigation, image fusion is used to provide additional information to the clinician in order to improve intervention outcomes [START_REF] Grimson | An automatic registration method for frameless stereotaxy, image guided surgery, and enhanced reality visualization[END_REF][START_REF] Cheung | Fused video and ultrasound images for minimally invasive partial nephrectomy: a phantom study[END_REF][START_REF] Shekhar | Live augmented reality: a new visualization method for laparoscopic surgery using continuous volumetric computed tomography[END_REF]. In this context, the image data is often acquired by different modalities and one of the input image pair suffers from poor signal-to-noise ratio. Moreover, registration time becomes an important aspect in a computer-aided navigation context.

In this work we focus on image registration in the context of partial hepatectomy which remains the preferable option to liver pathologies such as hepatocellular carcinoma, the primary hepatic malignancy with more than 780,000 deaths per year according to World Health Organization. During the hepatectomy, up to 50% of liver parenchyma containing the cancerous tissue can be removed. The resection must respect so called Couinaud regions which are 8 functionally independent segments, each having its own vascular inflow and outflow [START_REF] Strasberg | Nomenclature of hepatic anatomy and resections: a review of the brisbane 2000 system[END_REF]. Thus, accurate localization of both tumors but also vascular structures is necessary to perform the liver re-section. This localization is performed pre-operatively, based on contract-enhanced CT or MRI abdominal images [START_REF] Meinzer | Computerized planning of liver surgery-an overview[END_REF]. Beside helping the surgeon to create a mental image of the organ structure and its pathology, the geometry extracted from the pre-operative data can be visualized during the operation using the augmented reality [START_REF] Nicolau | Augmented reality in laparoscopic surgical oncology[END_REF].

However, the initial configuration of the liver at the beginning of the intervention can differ significantly from the preoperative state: in the case of laparoscopic surgery, large deformations are introduced by the pneumoperitoneum [START_REF] Tsutsumi | Image-guided laparoscopic surgery in an open MRI operating theater[END_REF]. For some cases in both laparoscopic and open surgery, it is recommended to perform the right hepatectomy in a semi-prone or a left lateral decubitus position (see Fig. 1a) in order to facilitate the access to the pathology [START_REF] Ikeda | Pure laparoscopic hepatectomy in semiprone position for right hepatic major resection[END_REF][START_REF] Wakabayashi | Recommendations for laparoscopic liver resection: a report from the second international consensus conference held in morioka[END_REF]. In [START_REF] Mulholland | Operative Techniques in Surgery[END_REF], the patient is positioned in full lateral position during right posterior sector robotic resection. Consequently, the pre-operative plan must be adapted in order to compensate for the deformations [START_REF] Jarnagin | Toward integrated image guided liver surgery[END_REF].

For this purpose, intra-operative imaging modalities are employed at the beginning or during the intervention [START_REF] Uchida | Recent advances in 3d computed tomography techniques for simulation and navigation in hepatobiliary pancreatic surgery[END_REF] and the plan is updated according to the acquired data. For example, if a contrast-enhanced CT is performed directly in the OR, the actual position of the pathology and vascular structures is localized directly in the image. However, this is rarely the case and usually, the intra-operative modality does not provide the sufficient image quality. In [START_REF] Kenngott | Real-time image guidance in laparoscopic liver surgery: first clinical experience with a guidance system based on intraoperative ct imaging[END_REF], a cone-beam computer tomography (CBCT) is used to perform intra-operative scan. CBCT is typically used in interventional radiology and its employment in during surgery remains rare. Nonetheless, this is changing thanks to the concept of a hybrid operating rooms which aims at proposing novel protocols that would benefit from the availability of various imaging modalities directly during the operation [START_REF] Benckert | The surgeon's contribution to image-guided oncology[END_REF]. While CBCT allows for reconstruction of the volume, the amount of information that can be extracted directly from the image is limited: in general, the liver tumors remain invisible under CBCT without intravenous contrast enhancement [START_REF] Chan | Target localization of 3d versus 4d cone beam computed tomography in lipiodol-guided stereotactic radiotherapy of hepatocellular carcinomas[END_REF]. However, CBCT is very sensitive to the acquisition timing [START_REF] Jones | Optimized dynamic contrast-enhanced cone-beam CT for target visualization during liver SBRT[END_REF] and according to [START_REF] Eccles | Intravenous contrast-enhanced cone beam computed tomography (IVCBCT) of intrahepatic tumors and vessels[END_REF], the intravenous injection does not guarantee the visibility of small tumors. While the internal structures are not easily visible under CBCT, it is possible to extract the liver surface using segmentation [START_REF] Li | Automated liver segmentation for cone beam ct dataset by probabilistic atlas construction[END_REF][START_REF] Li | Automated Liver Segmentation Method for CBCT Dataset by Combining Sparse Shape Composition and Probabilistic Atlas Construction[END_REF].

In this paper, we focus on the problematic of the intraoperative navigation during the partial hepatectomy performed in the lateral position to facilitate the access to the pathology. In this case, the resection plan based on the pre-operative CT scan acquired in the supine position must be considerably updated in order to take into account the large deformations of the tissue due to the supine-flank repositioning. We suppose that an imaging modality, such as CBCT, provides the information about the intra-operative shape of an organ, however, due to the low contrast and reduced X-ray dose, the actual locations of the internal structures and pathologies necessary to update the resection plan are not visible.

To this end, we propose a method that performs a fast reg-istration of the pre-operative data represented by a 3D model of liver and its internal structures and pathologies to the actual configuration given by the organ surface extracted from the intra-operative image. The method is based on a surfacematching algorithm which combines the iterative closest point technique with a biomechanical model based on a co-rotational formulation of linear elasticity in order to account for large deformations of the tissue. The algorithm thus requires extraction of sparse representations given by the pre-operative model and intra-operative surface point cloud from the input image data.

It is assumed that the sparse representation was acquired by a complete sampling, i. e. the intra-operative point cloud is distributed over the the surface of the organ.

The performance, robustness and accuracy of the method is assessed using a control and real datasets. In both cases, the target registration error (TRE) is computed inside the volume of the liver in order to assess the accuracy of the method which employs only the surface geometry to perform the registration. The control dataset is used to quantify the accuracy of the method via ground truth as well as to demonstrate the compatibility of the registration method with a real clinical scenario. For this purpose, the control dataset is represented by a semisynthetic pair of images, composed of a contrast-enhanced abdominal CT scan acquired in the supine position and a synthetic low-contrast CBCT volume mimicking the intra-operative flank configuration, generated using a simulation and image filtering. The extraction of the surface representation from the CBCT image is a time-critical procedure since its duration directly contributes to the time necessary to perform the intra-operative fusion. Therefore, this aspect is examined in details: two different methods of surface extraction are tested in terms of speed and registration accuracy to demonstrate the method can be performed in a real clinical setting.

The validation is further extended using nine real datasets, each composed of two 3D contrast-enhanced abdominal CT image acquired on a human subject scanned in lateral (flank) and supine positions. The goal of the real-data evaluation is to test the shape-matching algorithm in different patient-specific scenarios. Since in this case, the ground truth is not available, a novel validation platform is presented. It performs topological and geometrical matching of vascular trees, extracted semiautomatically from the supine and flank images, in order to construct a set of features distributed over the liver volume. The output of the validation pipeline is represented by the target registration error which is computed using the vascular features.

The proposed registration method is compared to a state-ofthe-art intensity-based technique employing free-form deformations [START_REF] Modat | A parallel-friendly normalized mutual information gradient for free-form registration[END_REF][START_REF] Modat | Fast free-form deformation using graphics processing units[END_REF] using the described validation platform. Then, both methods are examined in terms of parametrization; large-scale parametric studies are presented and it is demonstrated that a constant generic parametrisation can be found for the proposed method so that a close-to-optimal accuracy is achieved for all tested datasets.

State of the Art

In [START_REF] Crum | Non-rigid image registration: theory and practice[END_REF] and [START_REF] Rieder | Automatic alignment of pre-and post-interventional liver ct images for assessment of radiofrequency ablation[END_REF], registration methods are classified according to the type of the computed transformation as rigid, affine and deformable. The latter are employed in order to achieve an acceptable quality of fusion when the tissue undergoes large local or global deformations. In this case, the choice of the deformation model used by the method is important, as it defines the class of feasible transformations and impacts the number of degrees of freedom and parameters, as well as the computational time. In [START_REF] Maintz | A survey of medical image registration[END_REF], a classification of registration techniques is based on the type of information used in the optimization process. A large family of methods is based directly on voxel properties: they aim at maximization of a similarity metric computed using the image intensities, i. e. the dense representation of the input data. Another type of methods requires a sparse representation extracted from the image: the algorithms typically aim at minimizing an error metric given by the geometrical position of landmark features or the surface of an organ of interest. According to [START_REF] Maintz | A survey of medical image registration[END_REF], the methods using directly the voxel properties are considered as the most flexible, since they use all of the available information instead of relying on a sparse representation.

In the case of the intra-operative image fusion, the time needed to register the images becomes crucial as it is not possible to suspend the intervention in order to wait for the fusion results. Therefore, a method applied directly to the images should provide a more efficient tool. However, as we show in the paper, in the specific scenario with important tissue deformations, additional information represented by the biomechanical model and intra-operative surface of the organ may be necessary in order to perform reliable and fast fusion.

In what follows, we provide only a brief survey of non-rigid registration methods: since the number of publications reported yearly is huge, we do not aim at providing exhaustive survey of technique but we focus on those which are relevant w. r. t. the proposed method and its validation.

Intensity-based Registration

The registration methods based on voxel properties of the input images are usually implemented as iterative procedures optimizing a similarity metric such as sum of squared differences (SSD) or normalized mutual information (NMI). While the SSD metric can be applied only to images acquired using the same modality, the NMI metric assumes only probabilistic relationship between the image intensities, thus allowing for multimodal inputs [START_REF] Thévenaz | Optimization of mutual information for multiresolution image registration[END_REF][START_REF] Mattes | Pet-ct image registration in the chest using freeform deformations[END_REF][START_REF] Crum | Non-rigid image registration: theory and practice[END_REF]. These methods are applied directly to the input images which significantly simplifies the registration pipeline. Since the minimization process is ill-posed, it is usually regularized with additional terms such as elastic energy [START_REF] Marami | Model-based deformable registration of preoperative 3D to intraoperative low-resolution 3D and 2D sequences of MR images[END_REF] or diffusion [START_REF] Popuri | Fast FEM-based non-rigid registration[END_REF], numerically solved by a finite-element mesh covering the entire volume of the image.

Free-form deformations (FFD) [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast mr images[END_REF] are often used to implement the penalization term which is interpolated over a regular grid using B-splines. In [START_REF] Elhawary | Multimodality non-rigid image registration for planning, targeting and monitoring during ct-guided percutaneous liver tumor cryoablation[END_REF], [START_REF] Tang | MR-guided liver cancer surgery by nonrigid registration[END_REF] and [START_REF] Chen | Non-rigid mr-ct image registration for mr-guided liver cancer surgery[END_REF], multi-modal fusions of CT and MR images are presented. The methods start with a global rigid or affine alignment followed by a FFD-based registration compensating for local deformations. The performance of the metric optimization is improved in [START_REF] Song | Multi-modality liver image registration based on multilevel B-splines free-form deformation and L-BFGS optimal algorithm[END_REF] employing the L-BFGS method.

A modification of the mutual information (MI) metric is presented in de Manuel et al. ( 2014): the proposed organ-focused MI takes into account the mask of the liver obtained by segmentation. Recently, a hybrid method combining the metricminimization and landmark-based approach has been presented in [START_REF] Foruzan | Multimodality liver registration of openmr and ct scans[END_REF].

The main drawback of the registration techniques employing the FFD algorithm is the high computational cost. To address this issue, a parallel version is presented in [START_REF] Rohlfing | Nonrigid image registration in sharedmemory multiprocessor environments with application to brains, breasts, and bees[END_REF] running on a shared-memory supercomputer and in [START_REF] Rohrer | Parallel mutual information based 3D non-rigid registration on a multi-core platform[END_REF], a multi-core version is proposed. Finally, CUDA-based implementation is presented in [START_REF] Modat | A parallel-friendly normalized mutual information gradient for free-form registration[END_REF][START_REF] Modat | Fast free-form deformation using graphics processing units[END_REF]. It allows for choosing between two metrics (NMI and SSD) and provides four different regularization terms: bending energy [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast mr images[END_REF][START_REF] Rohlfing | Nonrigid image registration in sharedmemory multiprocessor environments with application to brains, breasts, and bees[END_REF], linear elastic energy, squared Euclidean distance of the displacement field and Jacobian determinant of the registration transformation [START_REF] Rueckert | Non-rigid registration using free-form deformations[END_REF]. The actual version is publicly available in a NiftyReg package, which is considered as the state-of-the-art implementation of the FFD registration. Although the method was originally designed for the deformable registration of the human brain, it was ranked in the top 5 algorithms in the Empire10 challenge dealing with the registration of pulmonary data [START_REF] Murphy | Evaluation of registration methods on thoracic ct: the empire10 challenge[END_REF] and recently, it was evaluated in an inter-patient registration of 20 abdominal data sets where it scored the best in terms of Dice metric [START_REF] Lee | Evaluation of Five Image Registration Tools for Abdominal CT: Pitfalls and Opportunities with Soft Anatomy[END_REF].

Geometry-based Registration

Unlike the intensity-based methods, the geometry-based techniques require additional preprocessing of the image data in order to extract the landmarks or surface of the organ. Since the optimization of the purely geometric error is also ill-posed, the process is regularized with additional a priori information which is often represented by a physical model of the tissue reconstructed from the extracted geometry.

A multi-organ deformable image registration used to analyze the motion of abdominal organs is presented in [START_REF] Brock | Accuracy of finite element model-based multi-organ deformable image registration[END_REF]. The algorithm minimizes the difference in surface shape and landmark displacement and it employs a mechanical model based on the finite element method [START_REF] Brock | Accuracy of finite element model-based multi-organ deformable image registration[END_REF]. In [START_REF] Vagvolgyi | Video to CT registration for image overlay on solid organs[END_REF], a shape-matching algorithm is used to register a pre-operative CT image to a point cloud extracted from stereoscopic intra-operative images. The method combines dense and sparse registration based on the iterative closest point (ICP) and point-based least-squares fitting, respectively. Only visual assessment of the algorithm is presented on kidney data.

A surface-matching approach was already employed in the context of intra-operative to pre-operative registration for hepatic interventions. In [START_REF] Cash | Compensating for intraoperative soft-tissue deformations using incomplete surface data and finite elements[END_REF], the ICP method is used together with a mixed co-rotational formulation of the linear elasticity in order to compensate for the tissue deformations: while the ICP is initially used to perform a rigid registration, the FE model is used to apply additional corrections compensating for elastic deformations. The reported fiducial registration error evaluated using a rubber-silicone liver phantom is 2 mm. A weighted ICP [START_REF] Clements | Robust surface registration using salient anatomical features for image-guided liver surgery: algorithm and validation[END_REF] is used in a rigid registration followed by elastic corrections to facilitate the tumor prediction and radio-frequency ablation [START_REF] Simpson | Model-Assisted Image-Guided Liver Surgery Using Sparse Intraoperative Data[END_REF]. The quantitative results confirming the benefits of the non-rigid registration are given only for phantom data. Similarly, liver phantoms are used to assess the role of the elastic compensation in non-rigid registration of pre-operative and intra-operative data in Rucker et al. (2013b,a). The wICP is further modified as the correspondences are updated at each iteration of the registration procedure. Also, the influence of the sparsity and partiality of the target data is studied w. r. t. the accuracy of the registration method. Finally, different acquisition modalities are compared in [START_REF] Wu | Registration of liver images to minimally invasive intraoperative surface and subsurface data[END_REF].

Recently, deformation correction algorithms have been evaluated on six patients undergoing open hepatectomy [START_REF] Clements | Evaluation of model-based deformation correction in image-guided liver surgery via tracked intraoperative ultrasound[END_REF]. Intra-operatively, tracked ultrasound probe was used to digitize the subsurface targets and laser range scanner was employed to partially reconstruct the surface. The authors report significant improvement of target registration error thanks to the deformation correction based on an elastic model reconstructed from the pre-operative data.

Deformable registration of intra-operative images acquired by a rotational C-arm during laparoscopic surgery to preoperative CT data is studied in [START_REF] Bano | Registration of preoperative liver model for laparoscopic surgery from intraoperative 3D acquisition[END_REF]. It combines rigid registration driven by anatomical features with deformable registration employing three landmarks. The accuracy of the method is evaluated using porcine liver data: the Euclidean distances between selected vascular bifurcations is reduced from 17.5±9.0 mm (max. 37.4 mm) before the registration to 10.3±2.7 mm (max. 15.8 mm). The method is extended in [START_REF] Oktay | Biomechanically driven registration of pre-to intraoperative 3D images for laparoscopic surgery[END_REF] by a specialized gas-insufflation model followed by diffeomorphic non-rigid registration and similar results are reported on synthetic data obtained by simulated deformation of porcine liver.

In [START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF], the non-rigid registration process is treated as an electrostatic-elastic problem. An electrically charged elastic FE model of liver [START_REF] Suwelack | Quadratic corotated finite elements for real-time soft tissue registration[END_REF] reconstructed from segmented CT data slides into a rigid target shape having the opposite charge. The target shape is reconstructed using the method proposed in [START_REF] Speidel | Intraoperative surface reconstruction and biomechani-cal modeling for soft tissue registration[END_REF] combining RANSAC calibration and correspondence analysis based on hybrid recursive matching. The method is evaluated using simulated data and a silicone phantom equipped with teflon markers. The initial mean Euclidean distance between markers is reduced from 23.9 mm (max. 46.6 mm) to 2.3 mm (max. 4.2 mm). The surface matching was also extended towards pose-independent scenarios [START_REF] Dos Santos | Pose-independent surface matching for intraoperative soft-tissue marker-less registration[END_REF].

A fast and reliable registration suitable for intra-operative scenarios where the tissues undergo large deformations remains a challenging task. Except for a few studies [START_REF] Vagvolgyi | Video to CT registration for image overlay on solid organs[END_REF][START_REF] Bano | Registration of preoperative liver model for laparoscopic surgery from intraoperative 3D acquisition[END_REF][START_REF] Clements | Evaluation of model-based deformation correction in image-guided liver surgery via tracked intraoperative ultrasound[END_REF], the assessment is performed or either synthetic or phantom data. Further, the robustness of the method w. r. t. low-quality intra-operative images is rarely discussed. Moreover, the parametrization of proposed methods often remains unclear. Nevertheless, if the method is to be applied directly during the operation, a patientspecific parametrization poses a serious issue to the usability of the registration procedure.

In the following section, we present the proposed registration method: Beside the detailed explanation of the iterative surface-matching algorithm, we also describe the data processing necessary to perform the fusion. The section 4 contains the validation of the method using the control semi-synthetic data. It describes the control data generation, presents the validation methodology and reports the results in terms of the initial and target registration error. The section 5 deals with the validation performed on the real datasets. The validation platform employing the vascular features is portrayed and the results for nine supine-flank pairs are presented for both the proposed method and an intensity-based registration method. In section 6, the parametrization of registration techniques is discussed. Then, parametric studies for the two methods evaluated on real data are presented together with consequences for the applicability of both methods in the clinical setting. Finally, the impact of the intra-operative data processing on the proposed registration method is briefly discussed and further perspectives are given.

Elastic Surface-matching Registration

Overview of the Registration Algorithm

The work flow of the surface-matching registration method is depicted in Fig. 2. The input data is represented by a pair of pre-and intra-operative images. While the former is a contrastenhanced 3D scan that allows for reliable localization of the surface and the internal structures of the organ of interest, the latter shows only the actual organ shape. Thus, due to a low contrast, the intra-operative image does not provide any reliable information about the position of the internal structures in the actual configuration.

The pre-operative data is used to reconstruct the biomechanical model: the liver parenchyma and the internal structures here represented by vessels are segmented from the input contrastenhanced images. This process can be done either manually or semi-automatically, but it must preserve a good level of detail of both the organ and the structures. Next, a finite element volume mesh M of the liver parenchyma is generated together with surface geometry of the vascularization. While the surface meshes are used only for visualization purposes after the registration is computed, the volume mesh M represents the domain of the reconstructed biomechanical model which is registered to the intra-operative data. In this paper, we always use the fully automatic method described in [START_REF] Boltcheva | Mesh generation from 3D multi-material images[END_REF] to convert a segmented map into a geometric structure such as a volume or surface mesh. The method first samples features in voxels located on the boundary of segmented region and then performs a Delaunay refinement to generate either the triangulation of its surface or a tetrahedralization of its volume.

Intra-operatively, the method first necessitates the extraction of the organ surface from the intra-operative scan. More precisely, a cloud C of points located on the boundary of liver parenchyma is constructed so that the points are distributed over the entire surface of the organ. This phase also requires segmentation of the liver parenchyma but unlike in the case of the pre-operative image, the intra-operative segmentation process must be done quickly and thus manual processing of a high-resolution intensity image is not possible. We employ a semi-automatic method proposed by [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF] which is applied to a down-sampled intensity image in order to speed up the segmentation. The algorithm is based on active contours and employs the level set method. As in the case of the pre-operative data, the point cloud C is extracted from the segmentation using the automatic sampling and refinement method [START_REF] Boltcheva | Mesh generation from 3D multi-material images[END_REF].

The last step of the intra-operative processing is the execution of the surface-matching algorithm registering the pre-operative model, the geometry of which is given by the FE mesh M, to the intra-operative point cloud C. In what follows, a detailed description of the algorithm is given.

Elastic Registration as Constrained Optimization

The proposed method of surface registration is formulated as an optimization problem min q Π(q) Φ(q) = 0

(1)

where Π(q) denotes the elastic energy of the pre-operative model and Φ(q) denotes the geometric constraints which impose the surface matching between the pre-operative model and intra-operative point cloud.

The Lagrangian L of the equality-constraint problem 1 is defined as

L(q, λ) = Π(q) -λΦ(q) (2)
where λ are the Lagrange multipliers associated with the constraints Φ. The minimization of the Lagrangian 2 is done by solving ∇ q,λ L = 0. The differentiation w. r. t. q and λ yields a system of equations

∇ q Π(q) -λ∇ q Φ(q) = 0 Φ(q) = 0. ( 3
)
The elastic energy of a deformable body is defined as Π(q) = 1 2 V σ εdV where σ and ε are the stress and strain in each position of the body, respectively. In this work, we employ the co-rotational formulation of the linear elasticity to define the internal forces g(q) = ∇ q Π(q) [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF]. In linear elasticity, the internal forces in a position q are defined as

gLIN (q) = B σ = B Dε = B DBu(q) (4)
where u(q) = qq 0 is the displacement of the particle at the actual position q w. r. t. its original position q 0 , and B and D are strain-displacement and stress-strain matrices, respectively. Due to the linearization of the strain tensor ε, the linear elasticity is not capable of handling large deformations. This issue is addressed by the co-rotational approach where the definition of the internal forces at position q is modified using a rotational matrix R which encodes the rotation of the local frame w. r. t. its initial orientation [START_REF] Müller | Stable real-time deformations[END_REF]. Then,

g(q) = R B DBRu(q). (5)
Since the rotation matrix R depends on the actual position q, the internal forces g are non-linear in q.

The continuous formulation is discretized using finite element (FE) method: we suppose that the FE mesh is composed of tetrahedral P1 elements equipped with linear interpolation functions defined in n nodes of the mesh. The vector of 12 discretized internal forces g e is computed using the weak form of Eq. 5 for each element e:

g e = v e R e B e D e B e R e u e ( 6
)
where u e are the actual displacements of the 4 vertices of the element e. For definition of other quantities, see Appendix.

Due to the non-linearity of internal forces, it is necessary to compute the linearization of g usually denoted as the tangent stiffness matrix. In the co-rotational formulation, this quantity is approximated for each element e as a 12×12 matrix

K e = v e R e B e D e B e R e . (7)
The local vectors of forces g e and element tangent stiffness matrices K e are assembled to obtain the global vector g composed of 3n values, and the global tangent stiffness 3n × 3n matrix K, respectively. The constraints Φ i (q), where i runs over m points of the surface point cloud C, are defined using a technique similar to the iterative closest point [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF]. For each node p i from the target cloud C, the closest point with position r i located on the surface of the mesh M is found. Since r i is located within a triangle (A, B, C) of the source mesh M, its coordinates are expressed using the nodal positions q A , q B , q C : r i = β iA q A +β iB q B +β iC q C where β iA , β iB , β iC are the barycentric coordinates of r i within the triangle (A, B, C). Denoting by n i the outer surface normal computed at position r i , the constraint Φ i associated with the cloud point p i is given as

Φ i (p i ) = (p i -r i ) • n i = (p i - v∈{A,B,C} β iv q v ) • n i . (8)
In order to solve the problem 3, it is necessary to compute the Jacobian of the constraints ∇ q Φ(q). Given the definition 8, this quantity is given as a m × 3n Jacobian matrix J which is highly sparse: on each line i, only nine entries are non-zero. Each non-zero triple corresponding to the triangle vertex w ∈ A, B, C is given as

J i,w = -β v n i .
Given the discretization and linearization of both elastic forces and constraints presented above, the Lagrangian system 3 can be solved by an iterative Newton-Raphson method where in each step, following augmented system of m+3n equations is solved for the positional correction ∆q and Lagrange multipliers λ:

K J J 0 ∆q λ = -g Φ . (9)

Iterative Solution with Numerical Damping

The matrix K itself is singular as no essential boundary conditions are applied after the assembly of the FE equations. At the same time, the definition of J does not guarantee the regularity of the augmented system. Therefore, we introduce a numerical damping by adding a small positive constant τ on the diagonal; see the Hessian modification methods in [START_REF] Nocedal | Numerical Optimization[END_REF]. The regularized augmented system corresponding to the optimization problem is then given by

K + τI K J J τI S ∆q λ = -g Φ (10)
where I K and I S are, respectively, 3n×3n and m×m identity matrices. The system 10 is solved by the method based on Schur complement as described below.

First, given the actual nodal positions q, the vector of internal elastic forces g e and the tangent stiffness matrix K e is computed for each element (steps 3 to 5 of Algorithm 1). Second, the local quantities are then assembled into global vector of internal forces g and the co-rotational tangent stiffness matrix K which is further modified by adding the damping factor τ on the diagonal (steps 6 to 8). Next, a vector ∆q un of corrections of the unconstrained positions is computed by solving the system K τ ∆q un = -g (steps 9 and 10). The vector of unconstrained positions q un = q + ∆q un is computed and employed to evaluate the actual constraints Φ and the corresponding Jacobian J. This computation requires very simple collision detection where for each point p i of the cloud C, the closest surface triangle given by vertices A, B, C of the mesh M is found, the point p i is projected on a location on the triangle (denoted as r i ) and the barycentric coordinates of posr i w. r. t. A, B, C are computed. This procedure together with evaluation of constraints and the constraint Jacobian is performed on lines 11 to 17.

In step 18, both K τ and J are used to compute the Schur complement S of the augmented system 9 as S = τI S -JK -1 τ J . In step 19, the dense m × m matrix S is used to obtain the vector λ of Lagrange multipliers; this is done by a projected Gauss-Seidel method solving the system Sλ = Φ. Finally, the vector λ is used to compute the corrected position update ∆q (step 20) which is applied to improve the actual estimation of the solution q (step 21 of Algorithm 1). The iterative process is stopped as soon as the the correction ∆q does not improve significantly, i. e., ||∆q|| < ε for a small constant ε.

Algorithm 1: Iterative Surface-matching Registration Input : Initial positions q 0 of source mesh M Positions p of points in target cloud C Output : Registered positions q f of source mesh M Parameters:

E [kPa], ν, τ 1 q ←-q 0 2 repeat 3 foreach e ∈ M do 4 g e ←-v e R e B e D e B e R e u e 5 K e ←-v e R e B e D e B e R e 6 g ←-e {g e (q)} 7 K ←-e {K e (q)} 8 K τ ←-K + τI K 9 ∆q un ←--K -1 τ g q un ←-q + ∆q un foreach p i ∈ C do r i ←-min x∈∂M ||x -p i || n i ←-normal of M in r i {A, B, C} ←-vertices of triangle containing r i {β A , β B , β C } ←-barycentric coordinates of r i Φ i ←-(p i - w∈{A,B,C} β iw q w ) • n i J i,w = β w n i for w ∈ A, B, C S ←-τI S -JK -1 τ J λ ←-solve(S, Φ) ∆q ←-K -1 τ Jλ q = q un + ∆q until ||∆q|| < ε q f ←-q 3.4

. Physical Interpretation of Numerical Damping

In this section, we analyze the impact of the numerical damping on the solution of the original problem given by 1. The damping is added to both the upper left part of the augmented system (via the term τI K ) and the lower right part (via the term τI S ). The former modification of the original system can be regarded as adding a lumped mass to elastic body. Thus, this modification does not influence the final solution of the system, but improves the stability of the iterative solution by increasing the inertia of the object being deformed.

The modification introduced by term τI S is more subtle. In case when τ = 0, the Lagrange multipliers λ, computed by resolution of system (-JK -1 τ J )λ = Φ, represent the forces necessary to impose the equality constraints 8, so that the surface of the deformed FE mesh matches strictly the points from the surface cloud C. The Schur complement JK -1 τ J has a physical interpretation: It represents the compliance in points r i of the biomechanical model located on the surface of the mesh, computed step 12 of the algorithm). The off-diagonal non-zero elements of the compliance matrix represents the mechanical coupling between the points r i , as they are attached to the FE mesh M. For more details about the compliance, see [START_REF] Peterlik | Constraintbased haptic rendering of multirate compliant mechanisms[END_REF].

In case when the damping factor τ > 0, the Schur complement is modified by matrix τI S . This modification can be also interpreted physically: It corresponds to the attribution of a non-zero compliance to each node p i of the point cloud C. Although this artificial compliance is not justified directly by the physical reality, as the point cloud C is a purely geometrical structure, it can account for inaccuracies associated with the extraction of the points in C. In other words, the non-zero compliance introduced by τI S results in softened matching between the surface of the FE mesh and the point cloud C.

Given this analysis, the damping coefficient τ can be regarded as a parameter complementary to the Young's modulus E: the values of E and τ influence the importance of the biomechanical model w. r. t. the constraints: for low values of both E and τ, the constraints are strictly respected, so that model is used mainly as an interpolator. If both E and τ are set to very high values, the stiff model ignores highly compliant target points resulting in very bad matching between the mesh and the surface cloud.

Validation Using Control Data

In the scenario considered in this paper, the intra-operative image represented by a non-injected CBCT scan is used to extract the surface point cloud necessary to perform the fusion with the pre-operative model. However, such data does not allow for quantitative validation because the non-injected intraoperative image does not contain any reliable information about the actual positions of the internal structures. Therefore, we perform two types of validation: in the first one, described in this section, we generate a synthetic CBCT image of a flank configuration from a contrast-enhanced abdominal CT scan of a real patient in supine position. The two images are then used as control data allowing for quantitative validation as described in this section. The second validation, presented in the following section, employs nine datasets of patients where both supine and flank positions are available.

The motivation of the validation on control data is twofold: first, we demonstrate the feasibility of the proposed registration method in the clinical setting. For this purpose, we evaluate the time needed to both process the input image data and perform the registration. The second motivation is given by the fact that in the case of semi-synthetic scenario, a quantitative assessment of the registration method can be performed using ground truth, represented by the simulated flank configuration.

Generation of Semi-synthetic Control Data

The simulated flank position is generated from a real contrast-enhanced abdominal CT image acquired by Siemens Somatom Sensation 64 scanner. The image was acquired in a study approved by the University of British Columbia Clinical Research Ethics Board (H08-02798). The generation process, illustrated in Fig. 3, consists of following steps: S1 The supine image is rotated by 90 • along the cranio-caudal axis to mimic the supine to flank rotation. The liver is segmented manually using an interactive tool1 . Two finite element meshes M INIT 1 and M INIT 2 are generated from the segmented map using the Delaunay refinement method presented in [START_REF] Boltcheva | Mesh generation from 3D multi-material images[END_REF] and available in a freelyaccessible library2 . The sizes of both meshes is given in Table 1.

S2 The fine tetrahedral mesh M INIT

1 is used to simulate the flank configuration. For this purpose, a gravitational force parallel to the coronal plane is applied. The boundary conditions are imposed so that they mimic anatomical constraints of the organ. The areas corresponding to the entry of the portal and hepatic veins are identified as well as the parts of the surface attached to ligaments. Since the latter cannot be extracted directly from the pre-operative CT, anatomical atlas was used to determine the location of ligaments approximately. Finally, the surface nodes located in the identified regions are fixed using homogeneous Dirichlet boundary conditions. The deformations are modeled using the MJED formulation of the hyperelastic StVenant-Kirchoff material [START_REF] Marchesseau | Fast porous visco-hyperelastic soft tissue model for surgery simulation: application to liver surgery[END_REF] available in an open-source simulation framework SOFA3 . The homogeneous model was parametrized with the Young's modulus E=5 kPa and Poisson's ratio ν=0.45. The coarse mesh is not used directly in the biomechanical simulation, nevertheless, it is deformed together with the fine mesh: at the beginning of the biomechanical simulation, each node of the coarse mesh M INIT 2 is mapped to an element of the fine mesh M INIT 1 using barycentric mapping. As M INIT 1 deforms by the simulation process, the coarse mesh follows its deformation because of the mapping; see for example [START_REF] Faure | SOFA: A multi-model framework for interactive physical simulation[END_REF] for further details. Thus, the simulation results in two deformed meshes M SIM 1 and M SIM 2 representing the same deformation. The initial and deformed configurations are depicted in Fig. 3b,c.

S3

The CT supine volume is warped using the displacement field provided by the difference of nodal positions of the initial fine mesh M INIT 1 and its deformed configuration M SIM 1 . We employ advanced warping method presented in [START_REF] Ručka | Fast reconstruction of image deformation field using radial basis function[END_REF] in order to obtain a realistic deformed image shown in Fig. 3d. Using the sparse displacement field given by the difference M SIM -M INIT , the warping algorithm first generates dense displacement field for each voxel of the image using radial basis function which is employed as interpolation inside the liver volume and extrapolation outside. Then, the warped image is generated by applying the dense displacement field to the original supine CT volume.

S4

The warped image still contains a high level of detail since it is generated from the contrast-enhanced CT. In order to mimic the intra-operative scenario, its quality is degraded in two ways. First, we suppress the contrast of Figure 4: The target registration error visualized using box plot (a,c) and histograms (b,d). On the left, TRE computed after registration employing the fine mesh M 1 and cloud C 1 extracted from the manual segmentation. On the right, TRE evaluated after registration using the coarse mesh M 2 and cloud C 2 obtained from the fast semi-automatic registration. The box shows the mean with standard deviation, the whiskers depict the 25th and 75th percentiles and the red crosses represent the outliers between the 75th percentile and the maximum TRE. Note the different range of y-axes in histograms.

the liver vascular structures: first, using the segmentation map of the vessels in the warped image, we replace the regions corresponding to theses structures in the intensity image by a 3D patch extracted from the non-vascularized parenchyma. Second, we apply the method presented in [START_REF] Sharp | Gpu-based streaming architectures for fast cone-beam ct image reconstruction and demons deformable registration[END_REF]

Validation Methodology

The semi-synthetic flank volume is used to assess the accuracy and feasibility of the registration proposed in 3. The validation procedure proceeds as follows:

Model reconstruction: The meshes M INIT 1 and M INIT 2 generated in step S1 described above are used by the biomechanical model employed by the registration method. We recall that only the fine mesh M INIT 1 was used to generate the synthetic flank image, while the coarse mesh was not directly employed in the simulation (S2). This mesh is utilized in the validation in order to avoid a potential bias introduced by using the same mesh to generate the synthetic data as well as to perform the registration.

Surface extraction: For the sake of validation, two clouds are extracted from the CBCT image generated in step S4. First, similarly as in the case of the pre-operative model, the liver is carefully segmented from the synthesized CBCT image using interactive manual segmentation and the surface mesh is extracted from the segmented map using the method based on Delaunay refinement; this time only surface patches are extracted since the volume elements are not needed. This process results in a point cloud C 1 composed of 594 points distributed over the surface of the liver.

Since a lengthy manual segmentation would not be feasible in the intra-operative clinical setting, we generate another cloud C 2 to demonstrate that the method is robust w. r. t. the segmentation quality. For this purpose, the CBCT image is first downsampled to 256 × 256 × 60 voxels, i. e. the number of axial slices is radically reduced in order to accelerate further processing. Next, the liver is segmented from the low-resolution image using a semi-automatic active-contour method [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF]. Due to the noise, it is necessary to apply several manual corrections, however, the time needed to obtain the low-resolution segmented map of liver including the necessary corrections does not exceed 10 minutes. Finally, the fast Delaunay refinement method is applied to the low-resolution liver resulting in a cloud C 2 composed of 576 points generated in less than 20 seconds. is the position of the same node the groundtruth mesh M SIM m . In order to quantify the improvement brought by the registration method, we also compute the initial error using the unregistered meshes M INIT m :

IE n = ||x SIM n -x INIT n || (12)
where and x INIT n is the position of the node n the initial unregistered mesh M INIT m .

Results

As described in section 4.2, the validation on synthetic data was performed on two different meshes M 1 and M 2 and with two different clouds extracted by manual and fast semiautomatic registration. The statistics (mean, median and maximum) of both the initial and target registration errors computed over all the nodes of the registered meshes are presented in Table 2. First, the initial error confirms important deformation of the liver: for nodes fixed by the boundary conditions, the initial displacement computed by the biomechanical simulation remains zero. Other parts of the liver subjected to the gravity force are significantly displaced up to more than 5 cm.

The reported results show that the registration method is capable of registering the tissue undergoing such significant deformation: the mean and median of the TRE is reduced under 4.2 mm and maximum TRE decreases to values close to 1 cm. Comparing the results obtained using different meshes reveals that the sensitivity of the method w. r. t. the 3D mesh used for the biomechanical model is negligible. The sensitivity w. r. t. the registration accuracy is higher, however, it remains acceptable: the deterioration due to the fast segmentation of low-resolution data of the both median and mean TRE remains under 0.5 mm while the maximum TRE is increased by less than 1.3 mm.

For the sake of better understanding of the error structure, the histogram of TRE is plotted in Fig. 4 for two extreme cases: the registration of fine mesh M 1 to cloud C 1 (extracted from the manual segmentation) and coarse mesh M 2 to cloud C 2 (extracted from the fast semi-automatic segmentation). The histogram shows that 90% of nodes are registered with error under 6 mm and TRE exceeding 8 mm is observed in less than 4% of nodes.

Further, a study of spatial distribution of the error was performed: the visualization over the fine mesh is visualized in several views in Fig. 5. The images indicate that the highest error is concentrated in the areas with a significant uncertainty related to segmentation of liver parenchyma. Typically, such areas are located around the entry of the portal vein: in noninjected images, it is often difficult to distinguish between the liver parenchyma and the vessel. Similarly, the error is higher in the left lobe which is very thin and it is difficult to distinguish the liver parenchyma from the stomach wall. This again leads to higher inaccuracy in segmentation, which in turn yields higher registration error. However, we recall that the goal of the method is to predict the position of internal structures inside the liver parenchyma, focusing mainly on the right part resected during the right hepatectomy. The reported results confirm that the method is capable of providing reliable information of displacements and deformations of these internal areas.

Finally, the registration algorithm required from 15 to 20 iterations of the damped Newton-Raphson method. The time needed to perform one iteration depends on the resolution of the mesh: 5 s in the case of M 1 and about 1 s in the case of M 1 . Together with the initialization of the FE model, the time needed to perform the surface-matching algorithm did not exceeded 120 s and 25 s for M 1 and M 2 , respectively.

Validation Using Real Data

In the previous section, we provided the validation of the entire registration process, focusing mainly on the intra-operative pipeline composed on the extraction of the sparse representation (the point cloud) from the intra-operative images and the surface-based registration. We assessed the accuracy of the method w. r. t. the quality of the input data and demonstrated that the method maintains the accuracy even if applied to data reconstructed by fast registration techniques compatible with clinical conditions.

In this section, we focus mainly on the surface-matching algorithm which is the core of the registration pipeline. We evaluate its robustness and accuracy on nine datasets, each composed of a flank and supine pair of contrast-enhanced CT images acquired on the same patient. The contrast enhancement introduced by the intra-venous injection allows for the extraction of vessel trees in both volumes. We employ topological matching to construct correspondences between the two trees, which are then used to compute the target registration error of the proposed surface-based method. The accuracy of the method is compared to a state-of-the-art intensity-based registration technique.

Validation Data

Nine data sets used for the evaluation were acquired in a study approved by the University of British Columbia Clinical Research Ethics Board (H08-02798). The contrast-enhanced CT images were obtained according to a current diagnostic protocol and each patient was scanned in the flank and supine position. The data sets were collected using three different scanners: Siemens Somatom Sensation, GE Medical Systems LightSpeed VCT and CT Toshiba Acquillion. Slightly different resolutions were used during the acquisitions: the voxel size ranges between 0.6 and 0.9 mm in coronal and saggital directions, and between 1.3 and 5 mm in the axial direction. Before further processing, the supine image was rotated by 90 • around the cranio-caudal axis to compensate for the rigid component of the body motion due to the supine-flank repositioning. Then, the two volumes in supine-flank pair were aligned according to an anatomical landmark: the bifurcation of the left and right portal vein which is close to the liver center of mass and is clearly visible in all data sets. After the initial rigid alignment, the sparse structures needed by the registration method were extracted from both volumes in each dataset. First, semiautomatic segmentation [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF] followed by manual corrections was used to segment the liver from each volume. Then, method based on the point sampling and Delaunay refinement [START_REF] Boltcheva | Mesh generation from 3D multi-material images[END_REF] was employed to construct the volume tetrahedral mesh from the supine segmented liver and the surface point cloud from the flank segmented liver. Exclusively for the validation purposes, the portal or the hepatic vascular tree was segmented manually from each volume. To facilitate the segmentation, the contrast of vessels was improved by the anisotropic-diffusion filter [START_REF] Enquobahrie | Vessel enhancing diffusion filter[END_REF]. The pairs of segmented trees were used to extract the vascular features for which a topological matching was used to establish correspondences between the flank and supine data for each dataset. Before describing the validation pipeline in detail, we present the algorithms for extractions of features and topological matching.

Extraction of Features from Hepatic Vascularization

The manual extraction of features from medical data such as the one presented in [START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF] usually provides only a small number of features. In the case of tissues with important vascularization such as liver, a set of features can be extracted from vessel bifurcations [START_REF] Lange | Validation metrics for non-rigid registration of medical images containing vessel trees[END_REF]. Nevertheless, with decreasing diameter of vessels, the number of bifurcations which can be reliably identified along a branch in both input images diminishes, although the branch itself is still visible. In what follows we describe an algorithm that overcomes this limitation.

First, the algorithm performs skeletonization of the segmented tree, i. e. it converts the image into a topological tree composed of branches represented by centerlines of the vessels. The skeletonization presented in [START_REF] Antiga | Centerline computation and geometric analysis of branching tubular surfaces with application to blood vessel modeling[END_REF] necessitates the construction of surface mesh, typically extracted from centerlines via marching cubes. However, the resulting skeleton often suffers from incompleteness as the method is very sensitive to the quality of the segmentation: as soon as the branch becomes too thin, the construction of the skeleton is interrupted and ignores the remaining part of the tree. Therefore, we employ a method based on the iterative Dijkstra minimum cost spanning tree presented in [START_REF] Verscheure | Three-dimensional skeletonization and symbolic description in vascular imaging: preliminary results[END_REF] and extended in [START_REF] Plantefève | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF]. The method is applied directly to the binary map of the vascular structure which is converted into a weighted 3D graph where each voxel belonging to the segmented volume is represented by a weighted node. The method constructs the skeleton in two steps: first, the spanning tree is constructed iteratively starting from the root voxel. The edges between voxels are constructed recursively using a sorted heap: in each step, the head of the heap having the minimal weight is marked and all its unmarked neighbors are inserted into the heap. The sorting weight of a voxel is defined as 1 r b which is the shortest distance of the voxel from the boundary of the binary map. This quantity is precomputed before the construction starts by an image filter. For each visited voxel, its distance from the tree root d r is stored.

In the second phase of the algorithm, the centerlines are extracted recursively from the spanning tree:

• The 0-order path P 0 is constructed as the shortest path connecting the root and the voxel with the highest d r in the graph.

• Any n-order path P n for n > 0 is extracted in three steps. First voxel expansion is performed, so that for each voxel r of the path P n-1 , a set V r of all voxels accessible from r is constructed. Second, a voxel t ∈ V v having the maximum value d r is found. Finally, the new path is constructed as the shortest path from r to t, denoted as (r, t).

After the skeletonization, a Bézier curve is fitted to each segment. The parametric Bézier representation has two advantages: first, it performs smoothing of centerlines which is important since the graph-based skeletonization often results in zig-zag shapes. Second, it allows for sampling an arbitrary number of points along the branches represented by continuous curves.

The skeletonization and curve fitting result in a skeleton given as a tree composed of smooth segments which are of two types: a closed segment connects two bifurcation points while an open segment is attached only to one bifurcation point; the other end of the segment correspond to the termination point of given branch. A schematic illustration of the skeleton is given in Fig. 6.

Topological and Geometrical Matching of Skeletons

The skeleton composed of parametric curves representing the vessel centerlines is constructed for both the flank and supine configuration of each dataset. Although both trees represent the vascularization of the same patient, their topology may differ due to the errors in segmentation which can miss some branches suffering from insufficient vesselness. We first perform topological matching between supine and flank skeleton in each dataset: we perform a semi-automatic method similar to that proposed by [START_REF] Lange | Validation metrics for non-rigid registration of medical images containing vessel trees[END_REF]. The method systematically tests possible matching hypotheses between the bifurcations by searching for corresponding paths. If a topological inconsistency is found, the hypothesis is rejected. If multiple hypotheses are accepted, they are examined interactively. Since the number of bifurcations is rather small (typically up to 20), the topological matching is fast and in spurious cases, it can be done manually.

The topological matching provides correspondences between the bifurcation points, thus it results in a coarse matching between segments which constitute both trees. Since the segments are represented by parametric curves, i. e. functions of a parameter t ∈ 0, 1 , it is possible to construct a fine matching as follows.

• Let S(t) be a closed segment of the supine tree and F (t) the topologically matching closed segment in the flank tree. Since the two segments match topologically, [S(0), F (0)] and [S(1), F (1)] are valid feature correspondences where the features are the matching end points (bifurcations) of the segments. Moreover, an arbitrary number of feature correspondences [S(u), F (u)] can be constructed straightforwardly for any chosen values of u ∈ 0, 1 .

• Let S(t) be an open segment in the supine tree and F (t) its matching segment in the flank tree. It is not important if F is closed or open: as S(t) is open, only one topological correspondence exists between the two segments: [S(0), F (0)].

In this case, it is not possible to perform fine matching using the sampling of parameters, because the terminal point S(1) does not correspond to the terminal point F (1); e. g. the supine segment is much shorter due to worst contrast of the input supine image. In this case, we propose to employ Euclidean distances to establish the fine mapping. We first discretize both samples by short lines. Next, we compute the length of each line and thus we approximate the total length of each segment. Let us suppose that S is shorter. We project S on F as follows: based on the approximative distance measure given by the discretization, we find e ∈ 0, 1 such that dist(S(0), S(1)) dist(F (0), F (e)). We then replace F with a new shortened segment F which corresponds to the fraction of F for t ∈ 0, e . We rescale the parameter of F so that F (1) = F (e). Then, S and F are treated as closed segments and the fine matching is performed as described above.

The fine matching illustrated in Fig. 6 results in a set of feature Table 3: Geometric characteristics of the data used in the surface-based registration and number of features utilized to compute the initial and target registration errors for both the surface-and the intensity-based methods.

pairs {(s f , t f )} where s f is the position of a feature extracted from the vascularization in the supine volume and t f is the position of (approximately) the same feature in the flank volume.

Validation Methodology

For each supine-flank pair, the validation employs the algorithms of feature extraction and matching as follows:

R1 Both the sparse representations needed for the registration procedure as well as the set of matched feature pairs are extracted from the flank and supine image as described in sections 5.1, 5.2 and 5.3.

R2

The liver mesh reconstructed from the supine image is registered to the surface point cloud extracted from the flank image using the surface-based registration presented in section 3. At the beginning of the registration, the positions {s f } of features extracted from the supine data are mapped using barycentric mapping to the undeformed mesh. As the mesh deforms during the registration process, these positions are displaced to {s f } which in fact represent the predicted positions of the selected features in flank configuration.

R3 For each feature f , the target registration error TRE is computed as

TRE f = ||s f -t f || (13)
where s f is the position of the feature f computed by the registration and t f is the position of the same feature extracted from the flank image. In order to quantify the initial difference between the supine and flank data, the initial error is computed using the non-registered positions s f of features:

IE f = ||s f -t f ||. (14)
The geometric characteristics of the FE meshes and points clouds, together with the number of features used to compute the error metrics are presented in Table 3.

Beside the evaluation of the registration error, we also present a comparison to an intensity-based method [START_REF] Modat | A parallel-friendly normalized mutual information gradient for free-form registration[END_REF] (b) . Therefore, we evaluate the target registration error TRE Nifty using the Eq. 14.

Results

The results presented in Table 4 show both the initial and target registration errors (IE and TRE, respectively) obtained for all flank-supine pairs. The table also reports the error denoted as TRE Nifty achieved by the intensity-based registration. Further details about the statistics of IE, TRE and TRE Nifty are visualized in Fig. 8.

First, the initial error IE evaluated on features extracted from rigidly-aligned data prior to the deformable registration confirms important deformation of the liver due to the supine-flank repositioning. Except for dataset P9, the mean initial error always exceeds 1 cm. The most significant initial deformations are observed in datasets P2, P6, and P7 where the mean and maximum error are close to or over 2 cm and 4 cm, respectively. A screen shot of a typical registration scenario is presented in Fig. 7 showing the flank and supine liver configurations before and after the surface-based registration.

The statistics of the target registration error TRE confirms the assessment performed on the semi-synthetic control data presented in section 4. In all cases, the maximum registration error is reduced under 1 cm and the mean error decreases under 4.5 mm. The evaluation of the intensity-based registration reveals, that the method is also capable of reducing the target registration error significantly. However, only in dataset P8 with the least significant deformations the values comparable to those obtained by the surface-based registration are achieved. In seven datasets, the maximum target registration error remains over 14 mm. In order to better understand the distribution of the error, the IE and TRE histograms are plotted in Fig. 9. In the case of the surface-based method, the TRE in 75% of features remains under 5 mm. However, in case of the intensity-based method, the same percentile is close to 15 mm.

The performance of the surface-based registration depends on the number of iterations of the damped Newton-Raphson method. Using the geometries listed in Table 3, the number of iterations never exceeded 60 and the time needed to perform single iteration remained under 1.2 s on a PC equipped with Intel i7-3770 CPU running at 3.4 GHz. Therefore, given the sparse representations of the source and target image data, the time necessary to perform the surface-based algorithm is remains under 70 s.

In the case of the intensity-based method, the time needed to compute the transformation heavily depends on the choice of the regularization terms: when computed on CPU (Intel i7-3770), the time varies between 1 and 30 minutes. For some regularization terms, a GPU version exists reducing the computation time to 2 minutes.

Discussion

Parametrization of Registration Methods

The important aspect of any registration technique is its parametrization: if the method depends on a large number of parameters and/or it is necessary to find the optimal patientspecific parametrization to achieve an acceptable accuracy of the registration, its deployment in an intra-operative scenario is highly questionable. First, in order to determine the parameters, it is necessary to have an access to the cost function used to quantify the optimality of given parametrization. In this paper, the cost function is represented by the TRE computed either from the ground truth in the case of synthetic data or from the vascular features extracted from both the supine and flank images for the sake of validation. However, in a clinical scenario, the cost function is typically unavailable. Moreover, optimizing for parameters requires iterative execution of the registration procedure. Again, this is not desirable in the intraoperative scenario as it would dramatically increase the time needed for execution of the method. Therefore, to make a registration method compatible with the intra-operative setting, it is necessary to find a generic parametrization, i. e. a parametrization which is not patient-specific but yields an acceptable performance for any dataset. The results presented in terms of TRE in Table 4 were achieved for optimal parametrizations reported in Table 5. These values were found via parametric studies executed for each data set. In the following section, we exploit the results of these studies to find the generic parametrization for each method.

Parametrization of the Surface-based Method

The proposed surface-matching registration method requires three parameters related to the biomechanical model and the solution method: Young's modulus E, Poissons's ratio ν and damping parameter τ described in section 3.2. Since a homogeneous model is employed, the Young's modulus and Poisson's ratio are constant over the volume.

We recall that in the proposed scenario, the role of the biomechanical model is to penalize the elastic energy of the deformation. During the registration, the deformation of the model is driven by positional constraints of the surface nodes and the set TRE mean (diff.) [mm] TRE max (diff.) observed quantity is also represented by positional data, the displacement of the points where the registration error is evaluated. Problems of this type are known as displacement-zero traction problems and in general, the solution of such problems does not depend on the absolute value of the elastic modulus [START_REF] Wittek | On the unimportance of constitutive models in computing brain deformation for image-guided surgery[END_REF]. However, in the proposed method, the surface displacements are not prescribed strictly due to the damping coefficient τ. This quantity has its own physical interpretation and plays a complementary role w. r. t. the Young's modulus as explained in section 3.4. In order to explore the influence of all the parameters on the accuracy of the method, we have performed a sensitivity study consisting of about 2,700 evaluation of the TREfor different parametrizations of the method. The parametrizations were generated systematically for all combinations of 4 different values of ν ∈ {0.2, 0.3, 0.4, 0.45}, 26 different values of E chosen from the interval 100 Pa, 50 kPa and 26 different values of τ chosen from the interval 0.001, 2.0 . First, we observed that typically the best registration results are obtained for Poisson's ratio ν ∈ {0.4, 0.45} being close to ν = 0.4 9 representing the incompressibility. Moreover, the overall influence of ν on the resulting error is very small: the variations of TRE max evaluated for different values of ν does not exceed 1 mm in all data sets and the impact on the TRE mean is less than 0.2 mm.

The sensitivity of the registration method is depicted in Fig. 10 where the maximum target registration error TRE max is plotted as a 2D function of Young's modulus and Poisson ratio. The first two plots showing the sensitivity in datasets P2 and P3 are typical representatives since very similar valley evoking the reciprocal function was observed also for P1, P4, P5, P7 and P9. In the case of P6, the shape of the function is slightly different as the accuracy deteriorates for high values of the damping coefficient τ in combination with low Young's modulus. This can be explained by higher sensitivity of given data set w. r. t. τ: the left bottom quadrant of the plot for P6 is similar to the shape observed for P2. Finally, the dataset P8 shows a different pattern: In this case, error surface is flat except for low values of τ = 0.001 in combination with E>10 kPa. This dataset is particular: the initial difference between the flank and supine con- figuration is relatively small when compared to other datasets and it is the only dataset where the surface-and intensity-based method give comparable results. These observations suggest that in the case of P8, the role of the model is limited to a simple interpolator and in that case, the method is insensitive to damping, thus, it is not necessary to account for the inaccuracies related to the segmentation. The parameter sensitivity studies confirm the complementary role of Young's modulus and the damping coefficient. Despite this complementarity, performing a parameter reduction by replacing the pair (E, τ) with a single parameter is not straightforward since while E is a parameter of the mechanical model, τ is a coefficient introduced by the iterative solution method. Further, the generic parametrization w. r. t. the maximum target registration error over the datasets d ∈ {P1, . . . P9} is proposed as follows. Let ε d (E, τ) denotes the TRE max of the dataset d, being defined as a function of the parameters E and τ. Further, let ε opt d are the optimal values of TRE max found for each dataset, i. e. the values reported in Table 4.

Then, the generic parametrization E gen and τ gen is computed by minimization arg min

E,τ d |ε d (E, τ) -ε opt d | (15)
i. e., E gen and τ gen are chosen so that the difference w. r. t. the optimal parametrization found for each dataset is minimized. The generic parametrization computed over nine supineflank datasets is given by E=3 kPa and τ = 0.02. The target registration errors obtained with this parametrization is presented in Table 6 for each dataset together with the difference w. r. t. the error achieved with the optimal parametrization presented in Table 4. Since the generic parametrization was constructed w. r. t. the maximum target registration error TRE max , this quantity increased for every dataset. The most significant difference w. r. t. the optimal parametrization is observed in datasets P7 and P9 but also in these cases, it does not exceed 2 mm. The differences in TRE mean are much less important: only in patient P7 the difference attains 1 mm; in other cases it remains under 0.3 mm and in P4 and P5 the generic parametrization (con-structed w. r. t. TRE max ) even slightly improves the mean target registration error.

Based on the available patient data, it is possible to conclude that using the generic parametrization to perform the registration of each flank-supine pair does not deteriorate the accuracy of registration. This result indicates that when it comes to parametrization, the method is compatible with a clinical intraoperative scenario.

Parametrization of the Intensity-based Method

We performed a parameter sensitivity analysis study also for the intensity-based registration method. In this case, the space was large due to the higher number of parameters given by 5 regularization weights listed in section 5.4 and the possibility to choose from two similarity metrics. Moreover, the execution time needed for single registration varies significantly between 1 and 30 minutes. For these reasons, it was not possible to perform the systematic search comparable to the one performed for the surface-based method.

Instead, we first performed the registration for 3,124 different parametrization sampled sparsely from the parametric space: each parametrization was given by the combination of all parameters set to one of the values from {0.0, 0.0001, 0.001, 0.01, 0.1}. We observed that some combinations resulted in over-constraining regularization, i.e., the registration transformation remained closed to identity. Alternatively, some combinations contained redundant terms: for example, in case of non-zero weights of the elastic terms, nonzero bending energy term does not influence the result significantly since the bending is already penalized by the elastic energy term. Based on this coarse initial study, we have selected tree pairs of regularization parameters for which we performed systematic sweeps: (w e1 , w e2 ), (w B , w J ) and (w B ,w L ) and for each parameter pair we computed registration for combination of parameters. Elastic and Jacobian energy weights w e1 , w e2 and w J were sampled over the entire admissible interval 0, 1 with step of 0.02. Since the registration accuracy exhibited much higher sensitivity w. r. t. w B this parameter was sampled from the interval 0, 0.0015 with step 3.10 -5 . andw. r. t. to Jacobian and bending energy regularization parameters (w B , w J ) for datasets P2 (c) and P9 (d). The combination of parameters resulting in the optimal values of TRE max reported in table 4 is marked with +. In case of elastic regularization, the combinations for which w E1 + w E2 ≥ 1 simply maintain the unregistered form of the source image.

In all datasets, the highest accuracy was attained with normalized mutual information similarity metric. While the regularization with non-zero coefficients w e1 and w e2 of the elastic energy term gives better results in datasets P5, P6, P7, the regularization with non-zero w B and w J yields higher accuracy in P2, P3 and P4. For the rest (P1, P8, P9), comparable results were obtained for both types of regularization. The regularization employing the L 2 norm of displacement field was not capable to achieve the optimal accuracy in any dataset for any w L > 0.

The examples of TRE max visualized as a function of either (w e1 , w e2 ) or (w B , w J ) are visualized in Fig. 11. The plots indicate that the choice of optimal parametrization is highly patient-specific. Moreover, the irregular shape of the error surface does not allow for choosing a reliable generic parametrization yielding close-to-optimal accuracy in all datasets. When compared to the parametrization of the surface-based method where such generic parametrization was successfully established for nine supine-flank pairs, the necessity of a patientspecific parametrization of the intensity-based method reduces the compatibility of this technique with the intra-operative setting.

Surface-based Registration in Intra-operative Setting

Finally, we briefly discuss the deployment of the surfacebased method in a clinical scenario. Since the construction of the 3D mesh is done pre-operatively, the procedure is not time-critical. However, this is not the case for the surface point cloud which must be extracted rapidly during the intervention as the first phase of the registration pipeline. Given a segmented map of an organ, a high-quality cloud point can be extracted in less than 20 seconds using the Delaunay refinement algorithm. Similarly, the surface-matching algorithm applied to the biomechanical model and the surface cloud point takes about 1 minute. Therefore, the time needed for the intra-operative segmentation remains the only critical factor.

In section 4, we demonstrated that liver can be segmented semi-automatically from the simulated CBCT volume in less than 10 minutes using low-resolution image, while preserving the acceptable accuracy of the surface-based registration. We believe that this procedure could be further improved in terms of both the speed and accuracy. For example, a fully automatic atlas-based segmentation of CBCT images is proposed in [START_REF] Li | Automated liver segmentation for cone beam ct dataset by probabilistic atlas construction[END_REF][START_REF] Li | Automated Liver Segmentation Method for CBCT Dataset by Combining Sparse Shape Composition and Probabilistic Atlas Construction[END_REF]. Further, the image acquisition systems including the intra-operative scanners make use of algorithms of automatic segmentation [START_REF] Heizmann | Assessment of intraoperative liver deformation during hepatic resection: prospective clinical study[END_REF]. Finally, the recent advances in machine learning show promising results in the area of automatic segmentation of medical images [START_REF] Christ | Automatic liver and tumor segmentation of ct and mri volumes using cascaded fully convolutional neural networks[END_REF].

Conclusions and Future Work

In this paper, we presented a model-based registration method based on a combination of co-rotational formulation of linear elasticity and iterative closest point. The method was proposed in the context of surgical navigation where an intraoperative modality provides the information about the actual shape of the organ, but it gives no information about the location of the important internal structures. As the case study, we focused on right hepatectomy which requires an update of the resection plan taking into account the supine-flank repositioning of the patient which results in large deformations of the liver parenchyma.

After describing the registration method together with necessary data processing in section 3, we first provided a validation employing control data 4. The goal of the assessment was twofold: (i) to quantify the registration error using a ground truth and (ii) to demonstrate the compatibility of the proposed method with the clinical setting. First, we presented a detailed description of the control data generation where a real contrastenhanced CT image acquired in supine position is transformed to low-contrast CBCT image mimicking the flank position. Then, two different FE meshes having different size of elements were reconstructed from the pre-operative image and two different point clouds were extracted from the synthesized image using manual and fast semi-automatic segmentation. Finally, the registration was performed for four combinations of meshes and point clouds in order to evaluate the impact of the input data quality on the accuracy of the method. Since the ground truth represented by a deformed model used in the generation of the CBCT image was available, initial and target registration errors were computed in the nodes of the FE mesh employed by the biomechanical model. Despite very large deformations yielding the maximum and mean initial error exceeding 13 mm and 50 mm, respectively, the method was capable to perform a physically valid registration: The mean and maximum target registration error were reduced under 4.1 mm and 11.1 mm, respectively.

In the second validation procedure, the error of the registration method was evaluated on nine pairs of real supine-flank images. In this case, a dedicated validation platform was proposed, based on a topological and geometrical matching of vascular trees extracted from the input images. Similarly as in the case of control data, the initial and target registration errors were computed. It was observed that in all patients, the maximum and mean target registration error were reduced to less then 10 mm and 4.5 mm, respectively. The method was compared to an intensity-based technique which aims at minimizing a similarity metric. However, only in one supine-flank pair, the results comparable to the surface-matching algorithm and in datasets with large initial error (mean > 17 mm and max. > 30 mm), intensity-based method was not capable to perform a reliable registration with acceptable accuracy.

Finally, a discussion dealing with the parametrization of the registration method was presented. It was shown that while in the surface-based method, a generic parametrization could be found, yielding a close-to-optimal for all tested datasets. Similar result was not achieved for the intensity-based registration where the parametrization seems to be patient-specific, thus the applicability of the method to intra-operative scenarios remains questionable.

As future work, we will aim at shortening the time needed to extract the intra-operative surface. We believe that using the state-of-the-art machine learning techniques, the extraction of the surface cloud can be significantly accelerated.     

where λ and µ are Lamé coefficients defined in terms of Young's modulus E and the Poisson ratio ν as

λ = Eν (1 + ν)(1 -2ν) µ = E 2 + 2ν . ( 19
)
In the co-rotational approach, it is necessary to compute the rotation matrix r e of the element e given its actual deformation. This is done by first computing the transformation between the original shape of the element and the actual shape given by 3 edges of the element attached to the same node. This transformation contains both the rotation and the deformation of the element and the rotational component is obtained using SVD decomposition which yields the 3×3 rotation matrix r e which is then organized into 12×12 matrix R e which consists of 4 copies of r e along its diagonal. For further details see [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF].

Putting it all together, the element internal forces are given by a vector of length 12 computed as

g e = V e
In both definitions, the integral over the volume of the element is computed as multiplication by the element volume v e , since the integrand does not depend on the position within the P1 tetrahedron thanks to the linearity of the interpolation functions.

  Figure 1: (a) Illustration of the lateral (flank) position during right hepatectomy. (b) Pre-operative contrast-enhanced CT scan acquired in the supine position. (c) Intra-operative CBCT scan acquired in lateral position.

Figure 2 :

 2 Figure 2: Scheme of the proposed surface-based registration: the biomechanical model of the organ including its internal structures is reconstructed from the pre-operative contrast-enhanced image (right). The surface point cloud is extracted from low-contrast intra-operative image (left). The biomechanical model is registered to the point cloud using the proposed surface-matching method (middle).

Figure 3 :

 3 Figure 3: The procedure used to obtain the semi-synthetic data: The liver if extracted from the contrast-enhanced CT (a) by segmentation and mesh generation (S1), resulting in model (b). The biomechanical simulation (S2) is applied to deform the model (c). The red points illustrate the boundary conditions imitating the ligaments. The initial and deformed FE meshes are used to warp (S3) the original CT. Finally, the warped image (d) is subjected to filtering and CBCT emulation (S4), resulting in deformed low-contrast CBCT image (e).

Registration:Figure 5 :

 5 Figure 5: Four different views visualizing the spatial distribution of the target registration error over the liver volume.

Figure 6 :

 6 Figure 6: Left: the pipeline showing the extraction and matching of vascular features. The contrast-enhanced CT images (a) are used to extract the topological trees (b) via skelenotization and Bézier curve fitting which provide feature matching (c). Right: the scheme depicts the vessel tree matching: first, the correspondences between the bifurcations (in green) are established via coarse matching of tree topologies. Then, fine matching is performed for each close and open segment depicted in red and blue, respectively.

Figure 7 :

 7 Figure 7: Visualization of liver and vascular tree from dataset P2: (a) Initial configurations after rigid alignment. (b) Configurations after the surface-based registration

Figure 8 :Figure 9 :

 89 Figure 8: The statistics of the registration error computed in the points sampled from the vascularization. For each data set, the graph shows the initial feature displacements computed before the deformable registration (left red box), the target registration error computed after intensity-based registration (middle green box) and similarly the target registration error evaluated after the surface-based registration (right blue box). Each box shows the mean with standard deviation, the whiskers depict the 25th and 75th percentiles and red crosses represent the outliers between the 75th percentile and the maximum TRE.

Figure 10 :

 10 Figure 10: Examples of the sensitivity analysis of the shape-matching registration w. r. t. the Young's modulus E and damping coefficient τ for datasets P2 (a), P3 (b), P6 (c) and P8 (d).The combination of parameters resulting in the optimal values of TRE max reported in table 4 is marked with +. In all cases, the Poisson's ratio is set to ν = 0.45.

Figure 11 :

 11 Figure11: Examples of the sensitivity analysis of the intensity-based registration w. r. t. elastic regularization parameters (w E1 ,w E2 ) for datasets P1 (a) and P6 (b), and w. r. t. to Jacobian and bending energy regularization parameters (w B , w J ) for datasets P2 (c) and P9 (d). The combination of parameters resulting in the optimal values of TRE max reported in table 4 is marked with +. In case of elastic regularization, the combinations for which w E1 + w E2 ≥ 1 simply maintain the unregistered form of the source image.

R

  e B e D e B e R e u e dV = v e R e B e D e B e R e u e (20) and the element tangent stiffness is given by a 12×12 matrix K e = V e R e B e D e B e R e dV = v e R e B e D e B e R e .

Table 1 :

 1 : radio-graphs are digitally generated from the warped CT volume with suppressed vessel contrast and then, a 3D CBCT volume is reconstructed using Feldkamp, Davis and Kress algorithm[START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. The resulting CBCT image has a resolution of 512 × 512 × 240 voxels; an axial slice of the volume is shown in Fig.3e. Characteristics of the FE meshes and point clouds used in the semisynthetic validation.

	Beside providing the displacement field necessary to synthe-
	size the flank CBCT volume, the simulation performed in step
	S2 also generates the ground truth represented by deformed
	meshes M SIM 1	and M SIM 2	which are explored to assess the accu-
	racy of the registration as explained in the following section.
		mesh M 1 mesh M 2 cloud C 1 cloud C 2
	# points	8,373	1,266	594	576
	# tetrahedra	41,437	5,272	-	-

  The lower table presents the statistics of the target registration error computed for the combinations of meshes and clouds employed in the validation.

			initial error [mm]
	mesh	IE mean	IE med	IE max
	M 1	13.9	12.2	54.8
	M 2	14.4	13.0	52.4
		target registration error [mm]
	registration	TRE mean	TRE med	TRE max
	M 1 to C 1	3.3	3.2	9.8
	M 1 to C 2	3.8	3.7	11.0
	M 2 to C 1	3.4	3.3	9.8
	M 2 to C 2	4.1	4.0	11.1
	Table 2: Statistics of the initial error computed using the simulated (ground-
	truth) meshes M SIM m	and unregistered meshes M INIT m	in the upper table.

Table 4 :

 4 Statistics (in [mm]) of the initial and registration errors computed for each dataset as well as across all the datasets (bottom line). The results were obtained for an optimal parametrization presented in Tab. 5; the questions related to the parametrization of both methods is discussed in section 6

		initial	surface-based		intensity-based
	set IE mean IE max TRE mean TRE max TRE Nifty mean	TRE Nifty max
	P1	11.9	35.9	4.5	8.7		7.8	19.3
	P2	19.1	48.5	3.4	8.5		11.4	39.8
	P3	11.9	29.8	4.0	9.7		4.9	15.0
	P4	17.7	33.3	4.3	7.3		11.1	19.8
	P5	13.0	32.1	3.7	7.9		7.0	14.9
	P6	20.7	39.3	4.1	8.0		9.8	24.4
	P7	22.4	51.4	3.9	8.8		8.2	22.5
	P8	7.2	16.9	3.3	8.3		3.5	8.5
	P9	11.0	26.7	2.8	6.2		5.8	13.4
	P*	15.0	51.4	3.8	9.7		7.7	39.8
		Intensity-based method	Surface-based method
	set	non-zero weights		metric	ν	E [kPa] τ
	P1		w e1 = 0.25		SSD	0.45	4.0	0.03
	P2		w J = 0.075		SSD	0.45	0.5	0.09
	P3		w J = 0.75		NMI	0.30	1.0	0.09
	P4		w J = 0.125		NMI	0.45	9.5	0.06
	P5		w J = 0.6		NMI	0.45	2.0	0.02
	P6	w e1 = 0.05, w e2 = 0.8	NMI	0.45	8.0	0.01
	P7	w e1 = 0.05, w e2 = 0.25	SSD	0.40	8.0	0.03
	P8		w e1 = 0.15		NMI	0.45	0.1	0.001
	P9	w e1 = 0.05, w e2 = 0.85	NMI	0.45	3.0	0.1

Table 5 :

 5 The optimal parametrization for each supine-flank pair for both intensity based (left) and surface-based (right) methods. These parametrization have been used to achieve the results reported in Table4. For the intensity-based method, only the non-zero values of regularization weights are reported.

Table 6 :

 6 Target registration error computed by surface-matching method parametrized with the generic parametrization where for each dataset, E=3 kPa, τ = 0.02 and ν=0.45. The table also show the difference w. r. t. to the values of the registration error achieved for the optimal parametrizations presented in Table 4.

ITK-SNAP: http://www.itksnap.org

CGAL: http://www.cgal.org/

SOFA: http://www.sofa-framework.org

http://sourceforge.net/projects/niftyreg

Scientific Cloud LM2015085, provided under the programme "Projects of Large Research, Development, and Innovations Infrastructures". We also thank J.G.Guevara for his valuable help with the generation of CBCT images from CT volumes.

Acknowledgment

We thank the Natural Sciences and Engineering Research Council of Canada and the University of British Columbia. Computational resources employed for the parametric studies were provided by the CESNET LM2015042 and the CERIT

Appendix

Here, we show the definition of the vector g e of the element internal forces and the local matrix K e assembled for each element of the model mesh M in every step of the Newton-Raphson method. Given a tetrahedral P1 element e equipped with four linear interpolation functions ϕ v (q x , q y , q z ) = a v +