Amir Yavariabdi

Adrien Bartoli

Chafik Samir

Maxime Artigues

Michel Canis

Mapping and Characterizing Endometrial Implants by Registering 2D Transvaginal Ultrasound to 3D Pelvic Magnetic Resonance Images

Keywords: Endometriosis, Deformable Surface, Slice-to-Volume Registration, Fusion, Localization

come

Mapping and characterizing endometrial implants by registering 2D transvaginal ultrasound to 3D pelvic magnetic resonance images

Introdution

Endometriosis is a multifocal gynaecological disease in which the cells from the lining of the uterus appear and grow outside of its cavity, typically in the pelvic area [START_REF] Chamie | Findings of Pelvic Endometriosis at Transvaginal US, MR Imaging, and Laproscopy[END_REF]. The symptoms of endometriosis vary widely from patient to patient, with some experiencing no symptoms, and others experiencing the whole gamut. Endometriosis affects approximately 10% of women of reproductive age and may cause chronic pelvic pain, abnormal bleeding, dysmenorrhea, dyspareunia, and infertility [START_REF] Mangler | Value of Diagnostic Procedures in Rectovaginal Endometriosis[END_REF]. The most common site where endometriosis occurs is the ovary, followed in descending order of frequency by deep lesions of the pelvic sub-peritoneal space, the intestinal, and the urinary systems [START_REF] Grasso | Pelvic Endometriosis: A MR Pictorial Review[END_REF]. There is currently no cure for endometriosis, but there are treatments for pain and infertility related to endometriosis (medications, hormone therapy, and laparoscopic surgery).

Endometriosis is a progressive disease which frequently infiltrates other organs and tissues. Many patients with endometriosis experience an average delay of eight years Email addresses: amir.yavariabdi@gmail.com (Amir Yavariabdi), adrien.bartoli@gmail.com (Adrien Bartoli), chafik.samir@gmail.com (Chafik Samir), maxwell8lf@hotmail.com (Maxime Artigues), mcanis@chu-clermontferrand.fr (Michel Canis) from presenting with symptoms to diagnosis and treatment [START_REF] Bokor | Semi-invasive Diagnosis of Endometriosis[END_REF]. Note that diagnosis of endometriosis should not just reveal the existence of the disease, but should also include such information as depth of infiltration and exact location with respect to the patient's anatomy. Surgery benefits from this information, since completeness of excision highly depends on the precision of diagnosis. Even though findings at physical examination may be suggestive, an imaging test is mandatory for diagnosis and preoperative surgical planning. The imaging modalities that have been mostly used to determine the stage of endometriosis are TVUS, MR, and laparoscopy (see Figure 1). Unlike MR and TVUS, in laparoscopy, diagnosis is confirmed intraoperatively. Diagnostic laparoscopy is the most widely used staging system for subtle lesions. In the case of deep infiltrating endometriosis, the value of laparoscopy is limited, as it does not allow the surgeon to see the endometrial implants that are located inside or under the organs and to quantify the depth of infiltration. To cope with these limitations, it is important to construct an endometrial map preoperatively using both MR and TVUS. TVUS imaging (Figure 1 (a)) is a reliable method for detecting small endometrial tissues and estimating their depth of infiltration. TVUS images are obtained using a two-dimensional array transducer that locally deforms the organs' shape, due to the probe's pressure, and captures a planar image. The main TVUS limitations include their limited field of view and low signal to noise ratio. MR scanning (Figure 1 (b)) is used as a complimentary examination, since it reveals the patient's pelvic anatomy and large endometrial tissues, and gives high resolution and fine detail 3D images of the patient's pelvis. However, MR imaging cannot show the small endometrial implants and the depth of infiltration.

By comparing these two imaging modalities, the radiologists are faced with an overwhelming amount of information, but the visual interpretation of these images is not an easy task. It is a challenging task because: 1) there is a local deformation between the two modalities, 2) in many TVUS data, the implants' neighboring organs are only partially visible due to the small field of view, and 3) the tissues that are exhibited in both MR and TVUS images do not belong to the same slice. Therefore, to ease the task of the medical experts in both interpretation and decision making, we need to move towards more comprehensive visualization techniques. To achieve this, an automatic fusion between TVUS and MR images is needed to remove many of the hurdles involved in determining the best intra-operative plan and transferring it to surgery. This will reduce the trauma done to healthy tissues and to avoid under-cutting the implants, which may cause recurrence. The registration process must facilitate the transfer of the two aforementioned information from a TVUS image to an MR data. This augmented MR will then contain superior information on the pelvic anatomy, including the location and shape of endometrial implants and their depth of infiltration in the host tissue.

The registration and fusion between TVUS-MR data is a difficult technical problem due to the large disparity in gray-level intensities (multi-modality), different dimensionality, soft tissue deformations, and the limited field of view and low signal to noise ratio of TVUS images. Therefore, deformable feature-based registration is an important tool in this context. Before we describe our feature-based method to register these two modalities, we summarize related research from the current literature.

Related Work

Image registration techniques can be roughly divided into two categories: 1) intensity-based and 2) featurebased. A general survey of deformable image registration may be found in [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF][START_REF] Modersitzki | Flexible Algorithms for Image Registration[END_REF][START_REF] Goshtasby | Image Registration Principles, Tools and Methods[END_REF][START_REF] Sotiras | Deformable Medical Image Registration: A Survey[END_REF][START_REF] Rohr | Elastic Registration of Multimodal Medical Images: A Survey[END_REF]. A thorough categorization of deformable medical image registration techniques has been presented by Sotiras et al. [START_REF] Sotiras | Deformable Medical Image Registration: A Survey[END_REF].

Registering MR with Ultrasound (US) images is not an easy task as: 1) the available organs in both modalities have different shapes and 2) the initial correspondences between shapes are unknown. To tackle these difficulties, various 2D/2D or 3D/3D registration strategies have been proposed with respect to the US image dimension. When using a 2D ultrasound scanner, the state-of-the-art methods mostly assume that the 2D US slice is parallel to one of the MR slices. With this assumption, one is able to manually select the corresponding 2D MR slice from the MR volume. This can be achieved automatically with the use of an Electro Magnetic (EM) tracker attached to the US probe during the acquisition which will give the spatial position of the US slice. It is then possible to locate an approximate position of the US slice in the MR volume. However, EM trackers are not available in most hospitals in Europe, so that most of the methods in the state-ofthe-art require an expert to search for the best MR slice match to the corresponding US slice.

Intensity-Based Methods

Intensity-based US-MR registration is a tremendously challenging task and is rarely implemented. The main reason of this difficulty is that US contains a speckle image of tissue boundaries, whereas MR provides information on tissue density. Generally, intensity-based US-MR registration methods consist of minimizing an energy function of two terms: 1) a data term that consists of an intensity measure (e.g. Normalized Mutual Information (NMI) and Cross Correlation (CC)) and 2) a regularization term that evaluates the smoothness of the deformation field (e.g. diffusion, curvature, and elastic smoothness).

Among all intensity-based techniques, the methods relying on NMI have been widely used to align US images to MR data. For instance, Mitra et al. [START_REF] Mitra | Prostate Multimodality Image Registration Based On B-Splines and Quadrature Local Energy[END_REF] present a method to register Transrectal US with MR prostate 2D images which uses B-spline free-form deformations with a new procedure of computing the NMI. In their work, the NMI is computed from the texture images generated from the amplitude responses of the directional quadrature filter pairs. They show that the entropy between US and MR images is typically more than the entropy of texture images due to variations in the gray levels. However, their method is not applicable to our registration problem. This is because in contrast with their images which include the prostate and its internal structures with very similar patterns, ours include multiple organs, and each organ has completely different internal structures. Thus, in our case, applying the directional quadrature filter cannot provide homogeneous signal intensity.

The advantage of using an intensity-based method is that it can directly use the image intensity information, without segmentation or user interaction, and thus can be achieved fully automatically. This registration approach also leads to a high registration accuracy. However, this approach is very sensitive to large intensity differences between the two images. Figure 2 shows an example of an intensity-based registration using NMI and B-spline deformation [START_REF] Rueckert | Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images[END_REF] to align TVUS to MR images. It is clear that NMI-based TVUS-MR registration method fails as the joint probability density function is too sparse, and there are very few pixels to calculate mutual information robustly. Therefore, there is no enough functional relationship between the intensity values in both images.

Feature-Based Methods

To solve US-MR registration, researchers have mainly used feature-based methods. In this approach, points or other extracted geometric features are used to register images. This registration approach is versatile in the sense that it can be applied to any image, no matter what the object or subject is. Another advantage of this method is that since the set of identified points is sparse compared with the original image content, a feature-based method has very fast optimization procedures.

Reynier et al. [START_REF] Reynier | MRI/TRUS Data Fusion for Prostate Brachytherapy[END_REF] propose a surface-based registration method for prostate brachytherapy. In their method, the prostate from the MR and US volumes is first manually segmented and then, the segmented surfaces are used to manually establish point correspondences. Thereafter, the US surface is rigidly registered with the MR surface. Finally, elastic registration is used to estimate the deformation between the two modalities. The advantage of their method is its ability to model the deformation using an octree-spline [START_REF] Szeliski | Matching 3-D Anatomical Surfaces with Non-Rigid Deformations Using Octree-Splines[END_REF]. Optimization is done with the Levenberg-Marquardt algorithm. However, the required procedure of placing point correspondences to align the surfaces requires significant user interaction.

In order to decrease user interaction various semi-automatic and automatic registration methods have been proposed. Cosse et al. [START_REF] Cosse | Diffeomorphic Surface-Based Registration for MR-US Fusion in Prostate Brachytherapy[END_REF] propose a semi-automatic method to register US-MR prostate and rectum surfaces. The MR image is first segmented using graph-cut while the US image is manually segmented. Registration is then carried out in two steps. First, rigid ICP is performed using the surface of the rectum and the prostate. Finally, a de-formable demons algorithm by Pennec et al. [START_REF] Pennec | Understanding the Demons Algorithm: 3D Non-Rigid Registration by Gradient Descent[END_REF] is applied to distance maps resulting from globally registered surfaces. Mitra et al. [START_REF] Mitra | Multimodal prostate registration using thin-plate splines from automatic correspondences[END_REF][START_REF] Mitra | A Spline-Based Non-Linear Diffeomorphism for Multimodal Prostate Registration[END_REF] automatically register 2D US to 2D MR to assist prostate biopsy. The 2D MR slice from the MR volume which corresponds to the US slice is manually selected. They automatically segment the prostate contour and use a triangulation approach to generate point correspondences. Deformable registration is then solved with a Thin-Plate Spline (TPS). The drawback of this method is that it cannot handle concave shapes.

To register a 3D MR data to a 2D US image, methods in the state-of-the-art require the selection of a 2D MR slice and use 2D/2D deformable image registration. This is a limitation, since a US image generally-at least in our case-intersects with multiple MR slices as it is not parallel to any of the MR directions; the problem is thus clearly a 3D/2D registration.

3D/2D Registration

With respect to data dimensionality, generally, medical image registration is either 3D/2D or 3D/3D. The former registration is a vital task in medical applications. In a surgical and radiological context, the precise knowledge of position (e.g. position of the surgical instruments or malignant tissues) is very consequential. With the help of 3D/2D registration the experts can obtain such information, and benefit in easier and better guidance in surgery. Here, we provide a brief summary of some relevant methods in 3D/2D registration in the context of treatment, image-guided intervention, and diagnosis. Most of the existing methods try to register 3D pre-operative data such as Computed Tomography (CT), MR, CT Angiography (CTA), to the 2D intra-intervention data including projective X-ray, CT-fluoroscopy, and optical images. Although our problem context is different, some of these methods are relevant as they try to register a 3D image to a 2D one.

2D/3D registration may refer either to registration of projective data to 3D data (for example, a 2D projective X-ray image with 3D CT images [START_REF] Groher | Deformable 2D-3D Registration of Vascular Structures in a One View Scenario[END_REF][START_REF] Livyatan | Gradient-Based 2D/3D Rigid Registration of Fluoroscopic X-Ray to CT[END_REF][START_REF] Yao | Deformable 2D-3D Medical Image Registration Using a Statistical Model: Accuracy Factor Assessment[END_REF]) or to register a single tomographic slice to 3D volumetric data (for example, a 2D US slice with 3D CT images [START_REF] Estépar | Towards Real Time 2D to 3D Registration for Ultrasound-Guided Endoscopic and Laparoscopic Procedure[END_REF]). The first case, which is known as projective registration, tries to find correspondences between a 2D image and a projection of a volume to planes. Biesdorf et al. [START_REF] Szeliski | Matching 3-D Anatomical Surfaces with Non-Rigid Deformations Using Octree-Splines[END_REF] proposes a method for the quantification of the aortic arch which combines 3D model-based segmentation with intensity-based elastic image registration. Their approach is based on minimizing a functional corresponding to a segmented 3D cylindrical tubular model as a first term and a 2D registered circular cross section as a second term. A parametric spline model is used to ensure smoothness of the resulting displacement field. Finally, they assume the knowledge of the projection from 3D to 2D and used this information to circumvent a 3D volume. The main drawback of this work is that the parametric spline model makes the method lose its generality in function spaces. In our problem, the cutting plane for TVUS is unknown and the projection is not available. Furthermore, we have an unknown deformation between the 3D and 2D objects. We assume the availability of segmentation to estimate the cutting plane and the deformation between objects.

The second case, which is also known as slice-to-volume registration, seeks to find correspondences between a 2D image and a cross-section of volume by a plane or warped surface. This can be considered as an extreme case of 3D/3D registration where one of the images reduces to one slice. One of the approaches to solve slice-to-volume registration is to use a so-called compounding technique [START_REF] Barry | Three-Dimensional Freehand Ultrasound: Image Reconstruction and Volume Analysis[END_REF][START_REF] Coupe | 3D Freehand Ultrasound Reconstruction Based on Probe Trajectory[END_REF][START_REF] Rohling | 3D Freehand Ultrasound: Reconstruction and Spatial Compounding[END_REF]. Heldmann et al. [START_REF] Heldmann | A Variational Approach for Volume-to-Slice Registration[END_REF] try to register a CT volume with US slices using this technique. In this manner, the US slice is compounded into volume by interpolating and then a 3D/3D image registration method is employed. However, they practically show that registration using the compounding technique cannot provide reasonable results. They also reveal that matching a CT volume to artificially reconstructed volumetric US data does not provide comprehensive information for the surgeon. Therefore, various other approaches have been used to solve the slice-tovolume registration problem.

Helmann et al. [START_REF] Heldmann | A Variational Approach for Volume-to-Slice Registration[END_REF] propose a variational deformable slice-to-volume registration to align series of 2D MR slices with 3D MR volumes. They minimize an objective function made up from a distance term (sum of square differences) and a smoothing term (curvature-based) with respect to a 3D nonlinear deformation field by using Gauss-Newton optimization. In contrast to compounding techniques, they evaluate the distance of the images only on the two-dimensional manifold where the data is known. The vital task in their algorithm is the choice of the regularizer. They show that at least a second-order regularization is needed to successfully avoid kinks and to estimate smooth deformation. Fei et al. [START_REF] Fei | Slice-t-Volume Registration and Its Potential Application to Interventional MRI-Guided Radio-Frequency Thermal Ablation of Prostate Cancer[END_REF] investigate intensity-based registration methods to rigidly align real time interventional MR image slices with high resolution preoperative MR volume. Their algorithm uses for interventional MR-guided radio frequency thermal ablation of prostate cancer. They use multi-resolution and multi-start strategies with two intensity similarity measures to avoid local maxima. The multi-start strategy is used to restart the registration process with randomly perturbed parameters obtained from a uniform distribution about the initial transformation values at the current resolution being used. They employ MI and correlation coefficient as similarity measures. The correlation coefficient is used at two lower resolution whereas MI is used at full resolution. This is due to the fact that at low resolution, MI surfaces are noisy and contain many local maxima. However, at full resolution, MI has a sharper peak than correlation coefficient. Ferrante et al. [START_REF] Ferrante | Non-Rigid 2D-3D Medical Image Registration Using Markov Random Fields[END_REF] propose a deformable intensity-based mapping algorithm between a 2D MR slice and a 3D MR volume. They try to find a linear plane transformation and an in-plane deformation field. This is achieved by using Markov Random Fields.

All these methods rely on an intensity-based similarity measure. However, US-MR registration with an intensity similarity measure usually fails. This is mainly due to the fact that this approach assumes that the images sufficiently correlated which is not the case in our problem. Thus, to tackle this registration problem we have chosen to pursue and develop a method based on corresponding geometric features between two images. Dalvi et al. [START_REF] Dalvi | Fast Feature Based Multi Slice to Volume Registration Using Phase Congruency[END_REF] propose a feature-based slice-to-volume approach to rigidly register 2D MR images to 3D MR volumes of the human brain. Their algorithm extracts phase congruency image features that are then matched using classical ICP [START_REF] Besl | A method for Registering of 3-D shapes[END_REF]. However, in order to estimate the deformation between two modalities, establishing point correspondences just at each MR slice are not sufficient due to the large MR inter-slice spacing. To solve this problem, the boundary of organs in each MR slice can be segmented and then a 3D surface may be reconstructed. This fills the MR inter-slice spacing by geometric information. Therefore, we chose to use contour to surface correspondences.

Our Approach

The aim of this work is to localize and charectrize endometrial implants in MR data. Since the small implant is not visible in the MR data and both itself and its neighbouring organs can be deformed in the TVUS images, there is no automatic correspondences between common features in both modalities. Therefore, an expert manually extract the boundary of available soft tissue organs in TVUS and MR images. In a case of 3D MR volume, we reconstruct a set of 3D surfaces (see Figure 3 (a)), while for the TVUS images we obtain a set of curves (see Figure 3 (b)). Now our goal is to focus on the shape of objects-for instance, a surface from MR data and a planar curve from TVUS-and fuse the 2D TVUS image to the 3D MR volume based on these obtained structures. When comparing two objects in different dimensions, there are mainly two possible solutions. One is to project the higherdimensional object to the lower space and then compare with the other object. However, in our problem, the parameters of projection and the plane of intersection are The general framework of our deformable slice-to-volume registration method. The MR data exhibits the patient's pelvic anatomy, whereas the 2D TVUS image shows the small endometrial implants and its host tissues along with depth of infiltration. After obtaining the images, we must provide the required data to our algorithm which are: 1) segmented endometrial implant and its neighboring organs in 2D TVUS data and 2) 3D models of the corresponding patient's pelvic organs which are constructed from a set of parallel 2D MR segmented organs. Our method computes a non-linear transformation embedding the 2D TVUS image in the MR coordinate frame while establishing point correspondences automatically.

unknown. The other way is to extrapolate the lower dimensional object (e.g. plane curves) to the higher space (e.g. space curves) and accomplish a more comprehensive search. Even though this approach can be computationally expensive, we must chose this approach due to our available data and information.

The dissimilarity between shape of organs in the MR and TVUS images results from two different transformations. Firstly, there is a difference in patient position due to the fact that TVUS acquisition is performed in lithotomy position, while MR imaging is performed in supine position. This transformation, denoted by T , introduces a change in translation and rotation of the 3D MRI surface S. Secondly, there is a nonlinear deformation, denoted by φ, due to the TVUS transducer's pressure, resulting in the deformed surface φ (T (S)). Let γ be a curve observed using the intersection of a cutting plane with φ (T (S)). Formally, our problem can be defined as: Given γ and S, estimate a curve on the S that optimally corresponds to γ under deformations. Note that T , φ, and cross-section of the MR volume (see Figure 3 (c)) are unknown variables and we need to estimate them in the registration process.

We propose a deformable registration method that may be applied to register a 2D TVUS with a 3D MR volume. We find a cross-section of the MR volume by a smooth surface, representing the TVUS. The endometrial implants and their depth of infiltration can then be mapped into the reconstructed patient-specific organ model from the 3D MR volume. Our variational one-step deformable ICP method directly registers a set of 2D curves (from TVUS) to a set of corresponding 3D surfaces (from MR). It computes a dense deformation field embedding the TVUS do-main in the MR coordinate frame while establishing point correspondences automatically. From this registration, any information marked in the TVUS frame such as the boundary of soft tissue organs and endometrial implants may be embedded in the MR frame. This includes the TVUS image itself, and we can thus directly visualize in the MR the deformation induced by the probe at the time TVUS was acquired. Figure 3 shows a general framework of our deformable slice-to-volume registration. In this figure, the warped TVUS surface and the location of the endometrial implant are shown in the MR frame.

Our contribution in this work are twofold. First, we bring a methodological contribution via the idea of combining TVUS and MR images to both characterize endometrial implants and localize them accurately with respect to one patient's anatomy. This may improve diagnosis and surgery planning [START_REF] Giusti | Anatomical Localization of Deep Infiltrating Endometriosis: 3D MRI Reconstructions[END_REF]. Second, we bring a technical contribution as a novel one-step ICP derived in a variational framework, and handling multiple curves to surface correspondences, while estimating a deformable transformation. Our method markedly extends the current literature on ICP. Lastly, an important advantage of our method over point-based methods is that it is extremely mild in terms of operator time.

The remaining of this paper is organised as follows: section 2 proposes a new deformable feature-based sliceto-volume registration approach. Section 3 provides the results of our experiments and a discussion related to the quantitative and qualitative registration results and accuracy of the proposed method. Lastly, section 4 concludes and outlines future work.

Methodology

Material

The study is carried out on ten patients who have small endometrial implants in the pelvic area. For each patient, the procedure of collecting data is as follows: first, a TVUS image which includes both the endometrial implants and their neighbouring organs is selected for each patient by an expert. Then, 4 or 5 2D MR slices from the MR volume that closely correspond to the moving 2D TVUS image are chosen by the expert. The MR images have an average size of 400×400×5 with a voxel resolution of 0.5×0.5×5 mm 3 . We assume that the TVUS image resembles the middle MR slice of the MR volume. It is, therefore, reasonable to assume that the TVUS slice is parallel to the corresponding MR slice in the first iteration of our algorithm. The soft tissue organs (bladder, uterus, ovary, rectum) and the implant in the TVUS slice and in the MR slices are segmented by an expert. Then, we use the method proposed by Kels et al. [START_REF] Kels | Reconstruction of 3D Objects from 2D Cross-sections with the 4-point Subdivision Scheme Adapted to Sets[END_REF] to reconstruct a 3D surface from a set of 2D MR contours. An example is shown in Figure 4 (a) and (b).

Registration Procedure

We choose TVUS to be our moving image (M ) and MR to be our reference image (F ). Segmentation results in a set of curves γ i and a set of surfaces S i , representing the boundary of corresponding organs in both modalities. For ease of understanding and simplicity of derivation, we first assume that γ includes only one curve and S contains only one surface as in the example of Figure 3. Let Ω be a bounded open set of R 2 representing the TVUS plane, q an arbitrary point in γ ⊂ Ω, and φ ∈ W (Ω, R 3 ) a 3D dense deformation vector field representing the slice-to-volume registration by embedding the TVUS plane in 3D while deforming it. W is a Sobolev space of sufficiently smooth vector fields over Ω with appropriate boundary conditions.

The literature of 3D/2D registration shows that deformable two-step ICP is a reliable choice when the projection information are available (e.g. [START_REF] Groher | Deformable 2D-3D Registration of Vascular Structures in a One View Scenario[END_REF]). This is due to the fact that projection information can be extremely useful to decrease the search space of the problem. However as we will see in the experimental results, standard ICP fails in the case of slice-to-volume registration. Therefore, we propose a new ICP-based algorithm which reduces the dependence of the procedure on the initial estimate and improves convergence rate.

A Variational Formulation for Two-Step ICP

Feature-based slice-to-volume registration is a difficult problem, especially for non-linear deformations. Without prior knowledge the problem is ill-posed. As prior, we use a local curvature regularizer. By doing so, the resulting displacement field is kept smooth. In other words, we impose a bounded Laplacian vector field to ensure that the registration outside γ would be estimated by interpolation in a way that the overall transformation smoothly maps the moving image onto the reference image.

We propose a variational formulation of the two inner steps of ICP with curvature regularization:

step 1: closest point computation ζ(q) = argmin Q∈S d 2 (Q, φ (q)) (1)
This implicitly defines ζ ∈ C 1 (γ, S), a function that computes the closest-point ζ(q) on S.

step 2: deformation estimation

φ = argmin φ∈W λ γ d 2 (φ, ζ (q)) dq data term + (1 -λ) Ω φ 2 2 dX regularization term
(2) where is the Laplacian operator and λ ∈ [0, 1] is a smoothing parameter.

Both steps in ICP must minimize the error, and thus, ICP is guaranteed to converge to a local minimum [START_REF] Chui | A New Point Matching Algorithm for Non-Rigid Registration[END_REF][START_REF] Fitzgibbon | Robust Registration of 2D and 3D Point Sets[END_REF]. Two-step ICP mainly depends on the nearest-neighbour heuristic used in step 1 which establishes binary point correspondences. It is easy to see that the procedure of establishing correspondences in two-step ICP makes it more vulnerable to local minima [START_REF] Chui | A New Point Matching Algorithm for Non-Rigid Registration[END_REF]. The question that may immediately arise is how the two steps can be merged into a single step to improve the convergence rate [START_REF] Yavariabdi | Contour-Based TVUS-MR Image Registration for Mapping Small Endometrial Implants[END_REF]. To achieve this, we follow [START_REF] Fitzgibbon | Robust Registration of 2D and 3D Point Sets[END_REF] and use a distance transform.

A Variational Formulation in One-Step Using Distance Transform

We now present our formulation of the registration problem by combining the two steps of the classical formulation in a single step. Substituting equation (1) in equation ( 2) yields:

φ = argmin φ∈W λ γ d 2 φ, argmin Q∈S d 2 (Q, φ) dq + (1 -λ) Ω φ 2 2 dX (3) It is clear that d 2 φ, argmin Q∈S d 2 (Q, φ) = min Q∈S d 2 (Q, φ).
This can be interpreted as the fact that the cost function depends on the distance to the closest point but not on the closest point itself. This allows us to rewrite equation (3) as:

φ = argmin φ∈W λ γ min Q∈S d 2 (Q, φ) dq + (1 -λ) Ω φ 2 2 dX
(4) We observe that the data term now involves a Distance

Transform D ∈ C 2 (R 3 , R), since D • φ = min Q∈S d 2 (Q, φ)
by definition. The distance map algorithm essentially assigns to every voxel (grid point) the distance to the nearest point. This is extremely efficient, fast and simple. Figure 4 is an example of the distance map which is generated from a 3D MR surface. We can now combine D with equation (4), leading to:

φ = argmin φ∈W λ γ (D • φ) 2 dq + (1 -λ) Ω φ 2 2 dX E[φ] (5)
Equation ( 5) represents a variational problem with E as cost functional. The next step is to compute φ using calculus of variation.

Euler-Lagrange Equation and Numerical Approximation

A function φ that minimizes E must fulfill the Euler-Lagrange (EL) differential equation. Several numerical optimization algorithms [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF][START_REF] Modersitzki | Flexible Algorithms for Image Registration[END_REF] may then be applied to solve the resulting nonlinear equation. In equation [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF], the data term places constraints on the deformation field φ at the curve location. To make it more general and to permit its formulation as an EL equation we use an index function δ γ , with δ γ : Ω → {0, 1}, δ γ (X) = 1 if X ∈ γ and 0 otherwise. We rewrite the cost function in equation [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF] as:

E[φ] = Ω λδ γ (D • φ) 2 + (1 -λ) φ 2 2 dX (6)
Let φ 0 be an initial registration. We use an additive update rule:

φ k+1 = φ k + U (7)
where U ∈ W is an incremental 3D dense displacement vector field. We restate our problem as:

min U E φ k + U (8)
By substituting equation ( 6) into equation ( 8), we obtain:

E[φ k + U ] = Ω λδ γ D • φ k + U 2 + (1 -λ) U 2 2 + φ k 2 2 dX (9)
The distance transform D is nonlinear and can be approximated by its first order Taylor expansion around φ k :

E[φ k + U ] = Ω λδ γ D • φ k + ∇D • φ k U 2 + (1 -λ) U 2 2 + φ k 2 2 dX (10)
where

∇ = [ ∂ ∂x , ∂ ∂y , ∂ ∂z ].
A function U that minimizes equation ( 8) must fulfill its EL equation. This is written as a system of three fourth order elliptic Partial Differential Equations (PDEs) represented by a 2 × 2 symmetric matrix equation:

µ γ D • φ k + ∇D • φ k U ∇D • φ k dq + ∂ 4 U ∂X 4 = 0 (11
) where µ = λ 1-λ . To solve equation ( 11) numerically, we discretize the curve γ in N points q 1 , . . . , q N and Ω on the rectangular pixel grid. Note that δ q (X) = 1 if ∃j ∈ 1, • • • , N such that X = q j and 0 otherwise. We consider the unknown function U = [u 1 , u 2 , u 3 ] on a rectangular pixel grid. Therefore, the discretization of equation ( 11) leads to a system of three PDEs:

µδ q (X) D φ k (X) + D x φ k (X) u 1 (X) + D y φ k (X) u 2 (X) + D z φ k (X) u 3 (X) D φ k (X) + ∂ 4 U ∂x 4 + ∂ 4 U ∂y 4 = 0 for ∈ {x, y, z} (12)
where D is the derivative of D with respect to ∈ {x, y, z}.

To solve the PDEs, we use a finite difference scheme with second-order boundary condition U = 0. This leads to a large but sparse linear system that can be solved by Successive Over-Relaxation (SOR) [START_REF] Saad | Iterative Methods for Sparse Linear Systems: Second Edition[END_REF]. The advantages of iterative solvers like SOR are twofold. First, iterative solvers generally perform very well in discarding the higher frequency parts of the error within the first iterations. This behaviour is reflected in a good initial convergence rate [START_REF] Bruhn | Real-Time Optic Flow Computation with Variational Methods[END_REF]. Second, they are suitable for solving large linear equations [START_REF] Jamil | A Comparison of Direct and Indirect Solvers for Linear Systems of Equations[END_REF]. The displacement vector field obtained from equation ( 12) is used to find the intersection of the MR volume by a warped TVUS surface and to localize an endometrial implant. To make our algorithm more stable, we use an a priori planar constraint. This is to avoid convergence to a spiky surface which would harm registration.

For that purpose, we automatically extract a plane from the TVUS point sets for the first few iterations. We choose the closest plane to the TVUS point set in the least squares sense [START_REF] Ahn | Least Squares Orthogonal Distance Fitting of Curves and Surfaces in Space[END_REF]. We project the displacement field on this plane. This choice puts constraints on the warped surface by flattening the TVUS free-form surface. Applying this strong constraint on the original TVUS surface improves the convergence of our algorithm. After a few iterations, we relax this constraint: when approaching the solution, the surface has only low frequency features and readily converges

q k ← q + φ 0
% The 2D points q embedded in the 3D space; φ 0 is an initialization of φ 7: Compute D and ∇D from Q % D and ∇D are used in equation ( 12) 8: while ∇U > ε do Compute U by solving equation [START_REF] Reynier | MRI/TRUS Data Fusion for Prostate Brachytherapy[END_REF] 10:

q k+1 ← q k + δ q k U
% Update the TVUS space curve 11:

φ k+1 ← φ k + U % Update the deformation field 12:
if k < ι then 13:

Π ← least-squares plane projection from q k+1

14:

q k+1 ← Π q k+1 % Flatten the TVUS space curve k ← k + 1 18: end while to the correct solution. The pseudo code of our method is given in Algorithm 1 where ε is a pre-specified convergence threshold set to ε = 10 -5 .

Handling Multiple Surfaces

We now extend our framework to handle multiple pairs of curves and surfaces. In this case, we have m curves γ 1 , . . . , γ m and a set of m corresponding 3D surfaces S 1 , . . . , S m . We use each surface S i to compute an Euclidean distance transform D i . By summing over all curve/surface pairs, equation [START_REF] Modersitzki | Flexible Algorithms for Image Registration[END_REF] becomes:

E[φ] = m i=1 Ω λδ γi D i • φ 2 + (1 -λ) φ 2 2 dX (13)
where δ γi (X) = 1 if X ∈ γ i and 0 otherwise. Using EL we obtain a system of three fourth order elliptic PDEs:

µ m i=1 γi D i • φ k + ∇D i • φ k U ∇D i • φ k dq + ∂ 4 U ∂X 4 = 0 (14)
Equation ( 14) is discretized by a finite difference scheme and leads to a sparse linear system that is solved using SOR, and eventually leading to an algorithm very similar to Algorithm 1.

Results and Evaluation

We demonstrate the accuracy of the proposed method by quantitative and qualitative tests on semi-synthetic and real world data sets. The true point correspondences are not available for real patient data sets. However, we assess the results with quantitative and qualitative evaluations on semi-synthetic data sets, generated to faithfully reproduce patient data features. We conduct more than a thousand experiments using semi-synthetic data with artificial deformations to test various aspects of the methods. The results are very informative. The experiments on semi-synthetic data emulate realistic acquisition situations and reveal the expected registration error. We calculate three different quantitative error measures to validate the results. To demonstrate the applicability in real applications, we also conduct tests on ten patient data. To show the generality of our method, we apply it to four different soft tissue organs: the bladder, uterus, rectum, and ovary.

Tests on Semi-Synthetic Data

Quantitative error measures. The first measure is the mean square distance between the N transformed moving points and the corresponding ground truth points (GT ). This measure is called the Mean Square Error (MSE), and is given by:

M SE = 1 N N j=1 q j -Q GT j 2 2 ( 15
)
where q j are the position of moving points after registration and Q GT j are the position of the corresponding points from ground truth.

Since the MSE does not evaluate the shape of the 2D curves embedded in the 3D space, we used another error measure which does not take the position into account, but only assesses the organ's shape. Being able to compare the curve after registration with the corresponding ground truth curve is essential to verify registration accuracy. To compare curves, we must first find an appropriate representation. To achieve this, the 2D space curves are simply parameterized by their arc-length and 20 points are sampled uniformly. Then, in order to evaluate the similarity between curve pairs, we compare the angle formed by adjacent points. This measure is known as the Shape Error (SE) [START_REF] Groher | Deformable 2D-3D Registration of Vascular Structures in a One View Scenario[END_REF] and is given by:

SE = 1 N -2 N -1 j=2 arccos (q j-1 -q j ) (q j+1 -q j ) q j-1 -q j q j+1 -q j -arccos Q GT j-1 -Q GT j Q GT j+1 -Q GT j Q GT j-1 -Q GT j Q GT j+1 -Q GT j (16)
The SE measure plays a significant role in identifying the similarity between the deformed moving curve and the ground truth curve. A low value of SE means that the curve shapes after registration are very similar even if misplaced.

To quantify the local registration error, we use the Target Registration Error (TRE). It is described as the mean square distance between M corresponding points not used in estimating the deformation. The main challenge in calculating TRE lies in finding corresponding target landmarks. The best target landmarks are those generated by a mechanical device which allows for an accurate positioning of the anatomy. However, this is not achievable for our study. The next best option is to choose anatomical landmarks in the patient's body and the centroid of an organ. We opt for this type of target points to estimate the local registration accuracy to validate our result. TRE is measured by:

T RE = 1 M M j=1 q T AR j -Q T GT j 2 2 (17)
where q T AR j are the position of target moving points after registration and Q T GT j are the position of the corresponding target points from ground truth. The target points used in our experiments are centroids and 10 boundary points which are not used to estimate deformation. A low TRE value shows good local registration accuracy. Generating semi-synthetic data. We create a test collection of 11 different 3D surfaces with 10 degrees of deformation and 11 different intersection curves, containing approximately 1200 data in total. To generate semi-synthetic data, the boundary of soft tissue organs such as the bladder, uterus, ovary, and rectum are segmented from five MR slices. Then, we use [START_REF] Kels | Reconstruction of 3D Objects from 2D Cross-sections with the 4-point Subdivision Scheme Adapted to Sets[END_REF] to reconstruct 3D surfaces from the 2D segmented contours. The reconstructed 3D surfaces are deformed using advection [START_REF] Bennett | Transport by Advection and Diffusion[END_REF]. In this method, every single point on the surface has a certain value of potential vorticity. The points are transported by their local displacement field on the surface. As a consequence, the surface can be deformed in complicated manners. The advection method is based on Lagrangian calculation. The advection equation is a hyperbolic PDE that governs the motion of the 3D surfac S as it is transported by a known displacement vector field U :

∂S ∂t + U • (∇S) = 0 ( 18
)
where ∂ ∂t is the derivative with respect to time. In our experiments, U is randomly generated and smoothed. The surface is deformed through increasingly larger degrees of deformation. The rate of deformation varies between 4% and 40%. Moreover, since the TVUS probe can move freely in any direction, it is also important to cut different views of the deformed surface. Therefore, various intersection curves are obtained at eleven different cutting angles in the range between 0 • and 20 • . These comprehensive simulations are used to validate our method and to test its robustness to organ deformation and probe orientation. An example of generated semi-synthetic data is shown in Figure 5. Results. In two-step ICP, the point correspondence problem is solved by the nearest neighbour heuristic which results in significantly poorer performance than the proposed method. To demonstrate the idea, we test both two-step ICP and our method on the same examples. Since we map the warped free-form surface to the planar domain in onestep ICP for 10 iterations, we use the same constraint for two-step ICP.

The qualitative results are shown in Figures 6 and7. In these two figures, we use a single 3D surface with 4% and 40% average deformation, respectively. We use two surfaces, one is the deformed version (Figures 6 and7 (b)) of the other (Figures 6 and7 (a)). The intersection curves are obtained at 20 • cutting angles. Then, we assume that these intersection curves are parallel to the middle intersection curve of the original 3D surface (see the red curves in Figures 6 and7 (a)). This assumption is made to keep the system similar to real world data. In these figures, the green curves are the ground truth solutions.

In Figure 6, we deform a 3D surface with 4% degrees of deformation. We then use both one-step and two-step ICP to find the best transformation to align the template set onto the ground truth set. The registration results are depicted in Figure 6 (c). At this deformation level, we can visually observe that both one-step and two-step ICP provide good results. This is visually demonstrated by Figure 6 (d) in which we show a 2D view (X-Y view) of the curves before and after registration. Moreover, to assess the behaviour of the method in a quantitative way, we use all the semi-synthetic data sets with 4% average deformation and calculate errors. The results are shown in rows 1 to 3 of Figure 8 (a). The MSE, SE, and TRE show that our method converges to a better solution at a low level of deformation. The two step ICP becomes less accurate as the initial solution is chosen further from the ground truth. However, by comparing the SE, it is easy to see that both methods have very similar shapes to their corresponding ground truth. Consequently, from these figures we can conclude that both methods provide good results at this deformation level.

To make the problem much more challenging, we use a 3D surface with 40% rate of deformation. The results are depicted in Figure 7. The two-step ICP shows very poor performance. The MSE, SE, and TRE in rows 7 to 9 of the Figure 8 show that our method still converges to a acceptable solution at 40% average deformation. The MSE, SE, and TRE indicate that the two-step ICP fails at this level of deformation.

In many cases, we observe that the warped intersection curves in one-step ICP are very close to the ground truth and that the method quickly converges to the optimum. In comparison, the warped intersection curves in two-step ICP gets trapped in local minima and mainly fail. The main reason of failure is large initial errors. In addition, failure may happen after the planar constraint is relaxed. Therefore, the warped contour may contain high frequency bumps which results in misplaced correspondences, and the registration process fails. Figure 7 shows an extreme example of such cases. Two-step ICP finds wrong correspondences and fails to retrieve the correct solution. We observe that the high frequency bumps have a strong influence on the matching process, causing the solution to get trapped in a local minimum. From Figures 8 (a) and (b), we can deduce that one-step ICP is stable and less sensitive to the initial solution. Moreover, the one-step ICP error gracefully degrades as the level of deformation and the cutting angle increase.

Tests on Real Patient Data

Data from real patients with small implants are used to validate our methodology of combining TVUS with MR and to compare our method against two-step ICP. In our experiments, small implants are seen in TVUS but cannot be identified from MR data. Therefore, 2D TVUS images are used to detect small endometrial implants and to segment the visible soft tissue organs whereas MR volumes are used to create a 3D geometric model of the corresponding organs. A TVUS image is manually selected for each patient by a clinician and approximately corresponds to the middle MR slice. Then, the soft tissue organs such as the bladder, uterus, ovary, and rectum as well as the implant in the TVUS slice and in the MR slices are segmented by an expert. Based on the available data and information, the initial solution φ 0 for both algorithms is chosen as follows: 1) the TVUS contour globally matches the middle MR slice and 2) the 2D version of the proposed method with a very small smoothing parameter (λ < 0.1) is used to locally register the TVUS contour to the middle MRI slice contour. To show how our method improves diagnosis, we use 10 different patients to localize small implants.

Figure 9 shows registration results for multiple surfaces. Figure 9 (a) shows the 2D TVUS slice. In this figure, the ovary, bladder, and endometrial tissue are segmented using expert knowledge. The soft tissue organs (ovary and bladder) and endometrial tissue can be identified by red and yellow, respectively. Figure 9 (b) shows MR slices with slice thickness of 5 mm. In this figure, the boundaries of the ovary and bladder are segmented from four MR slices. Then, we reconstruct the 3D surfaces from the 2D segmented contours. The reconstructed 3D ovary and bladder surfaces are shown in Figure 9 (b). We then use both one-step and two-step ICP to find the best transformation to align the 2D TVUS slice onto the 3D MR surface. The registration results for the proposed method and two-step ICP are depicted in Figure 9 (c) and (d), respectively. Figure 9 (e) illustrates the overlap images where the middle belongs to the rendered MR image and the top and bottom belong to the TVUS image. The resulting displacement field (Figure 9 (g) and (h)) is used to map endometrial tissue shape onto the MR surface, see Figure 9 (i) and (j). The overlap image can visually indicate that the bladder cannot be aligned by two-step ICP which leads to a low registration accuracy and inaccurate localization. On the other hand, based on our quantitative evaluations on the semi-synthetic data sets, we expect high registration accuracy for our method.

Our registration accuracy has also been evaluated by a pelvic radiologist. A precise registration algorithm must accurately align all corresponding organs and pelvic anatomical structures across modalities. The expert uses these features to visually assess our registration method. We provide the following data to the expert: MR slices, 3D patient-specific organ model(s), TVUS planar image, embedded TVUS surface, rendered MR image, and embedded TVUS surface with 3D MR surface. Examples of these are shown in Figure 9. We setup a registration quality score where 5 is excellent and 1 is bad. The expert evaluation is given in Figure 10 (a). The dash line means that the expert could not evaluate the results. Figure 10 (a) shows that our method mainly provides high registration accuracy, except for patient 9. From the expert's point of view, the registration results for patient 9 could not be evaluated since all organs (bladder and uterus) and pelvic anatomical structures (cervix, endometrium, and vaginal wall) in the TVUS surface were to be compared with 4 different MR slices (see Figure 10 (b)).

Discussion

The main advantages of TVUS-MR slice-to-volume registration over 2D/2D image registration are twofold. First, this approach removes the assumption that the TVUS slice is parallel to the corresponding MR slice. This assumption is not true as a TVUS image generally is not parallel to any of the MR directions. Therefore, it is important to discard this assumption, since the goal of endometriosis treatment is to resect the endometriosis without harming the surrounding healthy tissue. However, to achieve this, we need very accurate procedures in order to precisely localize endometriosis and enhance preoperative surgical planning. Slice-to-volume registration is then needed. As mentioned before, each slice of the MR volume can just depict the location of large tissues, and combination of all these slices can help the surgeons to visually reconstruct 3D structures of the patient's pelvic organs. However, such a reconstruction is prone to large error, since the pelvis has complex 3D structures. So a computer-aided reconstruction system must be used to accurately reconstruct 3D models of the patient's organs. This simply enhances the anatomy of this complex part of the body. Thus, the second advantage of slice-to-volume registration is that the surgeons do not need to deeply investigate each imaging modality, since all the necessary information needed to make an accurate preoperative surgical planning are provided as a new volumetric image which removes redundant information. This new image contains the patient-specific organ model from 3D MR volume as well as shape, size, location, depth of infiltration of endometrial implants from 2D TVUS images, and resection lines. This new image is promising to ease the task of surgeons to make better decisions to avoid under-or over-cutting during surgery. We evaluate the proposed method's performance in comparison with two-step deformable ICP. The quantitative error measurements on semi-synthetic data sets provide an assessment of our algorithm in a simulated clinical context. To demonstrate the applicability in real applications, we also conduct ten tests on patient data. Qualitative and quantitative tests on semi-synthetic data and clinical data sets clearly depict the accuracy of our method. Experimental error analysis show that our method remarkably reduce the error compare to two-step ICP. Experimentally, we see that the proposed algorithm has 45% higher registration accuracy on average than two-step ICP, from which we can conclude that our procedure is much less sensitive to increment of deformation and slant and can provide reliable results.

Conclusion and Future Work

We have described a variational one-step deformable ICP that directly registers a set of 2D curves in the 2D TVUS image to a set of corresponding 3D surfaces in the 3D MR volume. This method computes a dense deformation field embedding the TVUS domain in the MR coordinate frame while establishing point correspondences automatically. Our method uses Euclidean distance maps resulting from MR surfaces. It handles multiple pairs of curves and surfaces. Qualitative and quantitative tests on semi-synthetic data and clinical data sets show that our method remarkably reduce the error compare to two-step ICP.

The average execution time of the process is 153.82 ± 11.74 s with MATLAB code. We expect a speed up of the execution time of the model with optimized C++ coding or GPU programming and may thus use in real-time for multimodal fusion of female pelvic images. Besides the execution time, from a clinician's point of view, two fundamental questions have not been answered. The first one is how much registering a 2D TVUS with a 3D MR can improve the surgery planning. The second one is what the influence is of mapping an endometrial implant into a 3D patient-specific organ model on the surgery.

Finally, in order to improve surgical procedures to accurately resect the endometrial implants during laparoscopic surgery, we may enrich the surgeon's video data using Augmented Reality [START_REF] Collins | Realtime wide-baseline registration of the uterus in laparoscopic videos using multiple texture maps[END_REF]. This can be done with overlaying the reconstructed 3D patient-specific organ model including the implant's location, depth of infiltration, and resection lines in the laparoscope's video frame.

Figure 1 :

 1 Figure 1: The most commonly used imaging modalities for detecting endometriosis.

Figure 2 :

 2 Figure 2: Registration result for an NMI-based registration method[START_REF] Rueckert | Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images[END_REF]. In the overlap image (d), the top right and the bottom left are the deformed TVUS data and also the top left and the bottom right are the MR data. In (d), the red curve shows the uterus's shape after registration whereas the green curve shows the uterus's shape in the MR image. The large disparity between the green and red curves shows that even on this very simple case, NMI-based registration cannot cope with type of data.

Figure 3 :

 3 Figure3: The general framework of our deformable slice-to-volume registration method. The MR data exhibits the patient's pelvic anatomy, whereas the 2D TVUS image shows the small endometrial implants and its host tissues along with depth of infiltration. After obtaining the images, we must provide the required data to our algorithm which are: 1) segmented endometrial implant and its neighboring organs in 2D TVUS data and 2) 3D models of the corresponding patient's pelvic organs which are constructed from a set of parallel 2D MR segmented organs. Our method computes a non-linear transformation embedding the 2D TVUS image in the MR coordinate frame while establishing point correspondences automatically.

  (a) MR volume (b) 3D MR surface (c) Distance map

Figure 4 :

 4 Figure 4: An example of generating a distance map. (a) is the MR volume with segmented object (green curve). (b) is the 3D surface reconstructed from 4 parallel 2D MR curves. (c) is the distance map which is generated from the 3D MR surface. We use the distance map to solve ICP efficiently.

Figure 5 :

 5 Figure 5: An example of generated semi-synthetic data. In (a), the original surface is deformed through increasingly larger degrees of deformation. The rate of deformation are 4%, 20%, and 40%. (b) shows two different intersection curves obtained at 0 • and 20 • cutting angles.

One

  (a) Original (b) Deformed (c) Result (d) X-Y View

Figure 6 :

 6 Figure 6: Registration results for the proposed method and two-step ICP. The average deformation is 4%. (a): shows the reference 3D surface with the ground truth (green curve) and initial moving (red curve) curves. (b): shows the deformed 3D surface with the intersection curve which is obtained at a 20 • cutting angle (red curve). (c): shows registration results. (d): shows the 2D view of the curves before and after registration.

One

  (a) Original (b) Deformed (c) Result (d) X-Y View

Figure 7 :

 7 Figure 7: Registration results for the proposed method and twostep ICP. The average deformation is 40%. (a): shows the reference 3D surface with ground truth (green curve) and initial moving (red curve) curves. (b): shows the deformed 3D surface with 20 • cutting angles (red curve). (c): shows registration results. (d): shows the 2D view of the curves before and after registration. Best viewed in colour and close-up.

Figure 8 :

 8 Figure 8: Comparison of the proposed method with two-step ICP on semi-synthetic data sets with respect to cutting angle (slant) and deformation. (a): shows the results for 4%, 20%, and 40% deformations. (b): shows the results for 3 different intersection curves which are obtained at 0 • , 10 • and 20 • cutting angles.

Figure 9 :

 9 Figure 9: Example for two surfaces. The registration results for the proposed method and two-step ICP are shown in left and right columns, respectively. (a): shows a TVUS slice in which the bladder, ovary, and implant boundaries are depicted in red and yellow, respectively. (b): 3D ovary and bladder surfaces. (c) and (d): show warped free-form surfaces. (e) and (f) show the overlap images where the middle belongs to the rendered MR image and the top and bottom belong to the TVUS image. (i) and (j) depict endometrial implant after registration on the MR surface.

  (a) Evaluations (b) TVUS on MR frame

Figure 10 :

 10 Figure 10: Expert's evaluation of the proposed method. We setup a registration quality score where 5 is excellent, 1 is bad, and dash line means not rateable. (a): shows expert's evaluation result. (b): shows the TVUS surface embedded on the MR frame for patient data 9. Patient 9 can not be evaluated because the embedded 2D TVUS surface must be compared with 4 different MR slices.

  Algorithm 1 Pseudo-Code of the Proposed Algorithm 1: Input: A point cloud from MR (Q) & a set of points from TVUS (q) 2: Result: Deformation field (φ) % Grid of points h max × w max to represent φ

	3: Initialization:
	4: k ← 0	% Iteration counter
	5: ι ← 10	% Number of planarity constrained iterations
	6: