<figure><img src="image_58.png" /><figcaption>**Fig. 2.22:** Measurement of the detection error between a correct contour (green) and an arbitrary contour (red): (A) Hausdorff distance = length of s. (B) Average curve distance = average length of s. (C) Dice coefficient = $s/(s+t)$.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>**Fig. 2.23:** Illustration of results: green curve shows the manually segmented organ as the reference, red curve shows the result of semi-automatic algorithm and yellow line shows the Hausdorff distance. (A) Witness patient (1.17 mm/pixel). (B) Patient with endometriosis (0.49 mm/pixel). (C) Patient with prolapse (1.17 mm/pixel).</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>Fig. 2.24: Different exceptions of detection caused by unwanted effect in images: (A) As rectum may be folded somewhere, wrinkles appear in image as black rings. (B) The injected gel and intrinsic material may produce blur parts in image. (C) Intrinsic material show spots in rectum.</figcaption></figure>

<table><caption>Tab. 2.6: Summary table for results of the error measures between manual and automatic segmentation – the three sub-tables present results of each organ in different pathological cases.</caption><thead><tr><th>Patient (19)</th><th></th><th>Endometriosis (8)</th><th>Prolapse (8)</th><th>Witness (3)</th><th>Average</th></tr></thead><tbody><tr><td rowspan="3">Bladder</td><td>DC(%)</td><td>91.4</td><td>89.5</td><td>95.1</td><td>91.2</td></tr><tr><td>HD(mm)</td><td>9.46</td><td>10.87</td><td>5.74</td><td>9.46</td></tr><tr><td>ACD(mm)</td><td>1.35</td><td>2.23</td><td>0.90</td><td>1.65</td></tr><tr><td rowspan="3">Vagina</td><td>DC(%)</td><td>93.0</td><td>83.8</td><td>84.8</td><td>89.4</td></tr><tr><td>HD(mm)</td><td>6.14</td><td>11.34</td><td>3.02</td><td>7.84</td></tr><tr><td>ACD(mm)</td><td>0.89</td><td>1.80</td><td>0.67</td><td>1.24</td></tr><tr><td rowspan="3">Rectum</td><td>DC(%)</td><td>83.4</td><td>82.2</td><td>78.6</td><td>82.2</td></tr><tr><td>HD(mm)</td><td>31.12</td><td>25.91</td><td>33.59</td><td>29.32</td></tr><tr><td>ACD(mm)</td><td>5.04</td><td>4.77</td><td>5.66</td><td>5.02</td></tr></tbody></table>
### 2.3.3 Application: Pelvic System Mobility

Female pelvic disorders have a large social impact, the diagnosis of which relies on a key indication: pelvic mobility. The normal mobility is present in a healthy patient, meanwhile the hyper-mobility can be a sign of prolapse and the hypo-mobility for endometriosis. The evaluation of pelvic mobility is based on medical images. However, the latter does not provide directly precise values of these indicators. Moreover, suspension devices play an important role in pelvic organ function, but can neither be observed on medical images. Our objective is to propose an image-based analysis tool for the quantitative evaluation of pelvic mobility and the shear strain which has an impact on suspension devices. Hence using the approach described in the previous section we have developed a tool based on an efficient and semi-automatic motion tracking of multiple pelvic organs: the bladder, vagina and rectum presented in dynamic magnetic resonance imaging sequences. The computed displacement and shear strain field provide important information on the quality of suspension devices between organs, (i.e. on the condition of fascia tissues), to the opinions of medical doctors, which are interesting for a fine diagnosis in the clinical context, for example the early diagnosis of pelvic organ prolapse and the localization of possible lesion zones before surgery.

#### Tracking Algorithm

**Overview** A Lagrangian description is used to define the motions of the principal pelvic organs under deformation. Following this description in continuum mechanics, for any particle of the object, the change of its position over the time is defined by a function of transformation, from its initial configuration (or reference) to the current configuration. In other words, at any instant, the computation of physical quantities such as the displacement is always associated with the same particle (or material point), compared to its initial configuration.

To achieve the objective of motion tracking, firstly, we initialize the B-spline geometric model at instant $t_0$ (the first image of the sequence) as the reference configuration. Secondly, the model is updated for fitting each following image in the sequence to form the configuration at instant $t$. The process is recursive: the model of instant $t$ is deformed to fit the next image at instant $t + 1$, and it is also based on an energy minimization formulation. Once we obtain the models correspondent to all the instants, a displacement and shear field can be computed by using the analytical formulation of the models and following the Lagrangian description.
Meanwhile, we need to deal with the major concern about the plausibility of computed displacements. Normally, the description of the movement of material points remains in a fix landmark related to the initial configuration. So far the computation is based on the parametric coordinates, depending on the control points and their relative shift. But since these control points are not located on contours and they can freely drift along the contour. Thus the landmark for calculating the position of material points is mobile.

The idea is to fix some points on corners and to restrict the movement to perpendicular to the contour. For this purpose, the B-spline model is divided into piecewise cubic Bézier curves joined end to end. The ending point of one Bézier curve is attached with the starting point of another to form a connection point which is considered as a material point and where the curvature is higher than other-where. Here impose we an anatomical assumption that the salient points on the organ surface (or corners on the contour in a 2D plan) represent material points, which is based on the observation of medical images. Such material points, denoted as features in the remainder of this paper, are firstly tracked by the optical flow algorithm [LK+81], which consist in the beginning configuration for each recursive process of deformable registration from t to t + 1.

**Features definition** We apply a new adjusted model (composite Bézier curve) and the connection points of Bézier curves become the material points. At the initial step, the choice of these features is performed automatically based on detection of the higher curvatures. Suppose that the contour of an organ in the initial configuration (instant $t_0$) is defined by a parametric B-spline curve $B^0(u) = \sum_{i=0}^{n-1} N_{i,3}(u)p_i$. The curvature is related to the first and second derivative of the parametric function, which can be calculated by using the standard parametric coordinates as follows:

$$ \kappa(u) = \frac{\|B^{0'}(u) \times B^{0''}(u)\|}{\|B^{0'}(u)\|^3} \qquad (2.24) $$

where $B^{0'}(u) = \frac{dB^0}{du}(u)$, $B^{0''}(u) = \frac{d^2B^0}{du^2}(u)$. Hence, the curvilinear curvature can be computed directly based on the initial B-spline model. In an automatic way, the algorithm enables to select a number of features where the curvature is higher by using the following criteria:

* MaxNbFeatures and MinNbFeatures maximal and minimal number of features to be selected;
<figure><img src="image_61.png" /><figcaption>Fig. 2.25: (left-middle) Choice of the features along the B-spline contour. (right) Adjusted model at instant $t_0$. The figure is generated by our application.</figcaption></figure>

* MinDetectRadius minimum distance between two adjacent features (used in order that the features are well spaced);

* CurvThreshold threshold of the curvature (used to avoid the features on flat segments).

Figure 2.25 (left, middle) shows the selected features with a curvature map bounding the contour. Meanwhile, the features are resorted with respect to the parametrization of the B-spline. Let $[\mathbf{pf}_0, \mathbf{pf}_1, \dots, \mathbf{pf}_l]$ be the selected feature points, where $\mathbf{pf}_i = B^0(u_{k_i})$. $l+1$ piecewise Bézier curves can be created and joined at these points, who represent the same geometry as the initial B-spline and constitute the adjusted model. As defined in many standard texts, the Bézier curve is a special case of the B-spline who has four control points. In this case, for each piecewise Bézier curve, the two extremities are fixed by two adjacent features. Let us take the first Bézier segment between $\mathbf{pf}_0$ and $\mathbf{pf}_1$ as an example. This curve, denoted $S_0$, is defined by the Bernstein polynomial $M_{i,p}: u \in [0, 1] \to [0, 1]$ as:

$$ 
\begin{aligned}
M_{i,p}(u) &= \frac{p!}{i!(p-i)!} u^i (1-u)^{p-i} \\
S_0(u) &= \sum_{i=0}^{3} M_{i,3}(u) \mathbf{p}_i
\end{aligned} 
\qquad (2.25) 
$$

where $p_0 = \mathbf{pf}_0$ and $p_3 = \mathbf{pf}_1$. Then the two free control points $p_1$ and $p_2$ are determined such that each Bézier curve fits the original B-spline contour on the interval $[u_{k_0}, u_{k_1}]$. The curve fitting method is well known as a standard approach by solving a system of linear equations. Figure 2.25 (right) shows the adjusted Bézier model as the initial configuration, created from the original B-spline model.

**Features Tracking** The approach of optical flow proposed by Lucas et al. [LK+81] can be well adapted into our needs. The method is based on the assumption of the constant pixel intensities and their low local variation. Because of the local
<figure><img src="image_62.png" /><figcaption>Fig. 2.26: Tracking of the features using the optical flow method: window grid size = 30 × 30 pixel.</figcaption></figure>

evaluation, the method would fail if the deformation between images was large or the intensity in a region was uniform. However, in our case, the features to be tracked are located on the edge of organs where the variation of pixel intensities is large. Moreover, with the acquisition frequency of about one image per second, the deformation between two images (instant t and t + 1) is relatively small. All these aspects permit a good performance of the optical flow tracking. In a practical view, it is necessary to choose the appropriate window grid size for the computation which respects the propriety of locality and where the image pattern is enough representative to be tracked. Generally, the choice of size depends on the image, we have adjusted this parameter empirically during tests (see Figure 2.26 for illustration of results). Once the tracking of features is completed using the optical flow, the transformation of each feature is to be applied over its two adjacent free control points to form the intermediate model at instant t* (between t and t + 1). Let always $S_0$ be the first Bézier curve with four control points $[\mathbf{p}_0 = \mathbf{p}f_0, \mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3 = \mathbf{p}f_1]$. After applying the transformation, we obtain the following formula:

$$ S_0^*(u) = M_{0,3}(u)T_0(\mathbf{p}_0) + M_{1,3}(u)T_0(\mathbf{p}_1) + M_{2,3}(u)T_1(\mathbf{p}_2) + M_{3,3}(u)T_1(\mathbf{p}_3) \quad (2.26) $$

where $T_0(\mathbf{p}f_0) = \mathbf{p}f_0^*$ and $T_1(\mathbf{p}f_1) = \mathbf{p}f_1^*$ determined by the optical flow method (see Figure 2.27 for illustrations).
<figure><img src="image_63.png" /><figcaption>**Fig. 2.27:** Deformation of the model from instant *t* to *t* + 1: firstly, the features are deformed using optical flow tracking method, for example the transformations of point 0 and 1 are computed by the optical flow equation. Then, the same transformations are applied on their neighboring free control points to form the intermediate model (instant *t*^*). Finally, the free control points (for example the ones situated between 0 et 1) are adjusted in the direction perpendicular to the segment using virtual image correlation, and produce the model at instant *t* + 1.</figcaption></figure>

**Optimization** The last step consists in adjusting the free control points associated with each Bézier curve so that the model fits perfectly the organ in images. In order to minimize the drifting along the contour, the displacement (transformation) of these points is restricted to the direction perpendicular to the segment linking two adjacent features. And the magnitude of the displacement is also limited to be small because physically no new salient points will be produced during the deformation. Moreover, after the step of transformation based on the optical flow, the intermediate model has been already close to the deformed organ in the image. However, the adjustment allows a finer fit between the model and the image at instant *t* + 1. Following the same principle as the B-spline based detection, the adjustment of positions is carried out by minimizing the difference between the virtual image and the real one.

**Displacements and Shears** Our approach relies on a Lagrangian description which enables to compute continuous displacement field with our analytical geometric models. Following this description, the displacement field of each organ must be evaluated separately and each organ must have a reference configuration. This reference corresponds to the shape of the organ wall found at instant *t*₀ (on the first image, detected by our registration method). However, practically a given number of discrete points are sampled in the zone of interests in the initial configuration, and the analysis of mobility will rely on this discrete domain. The computation
<figure><img src="image_64.png" /><figcaption>Fig. 2.28: Computation of the displacements (A) and shears (B), between instant $t_0$ and $t$.</figcaption></figure>

of displacements and shears is then based on the parametric coordinates of Bézier curve

After computing the displacements for each Bézier curve respectively, we obtain a set of vectors $dp_i^t$ which defines the displacements field at instant $t$, and this computation is repeated for each instant (image) of the dynamic MRI sequence. Concerning the shear strain between organs, the middle organ (vagina) is chosen to be the reference organ. Analogous to the computation of displacements, at the initial configuration we generate several couples of discrete points from different organs based on the zone of interest defined by users. These couples of points are linked by segments perpendicular to the vaginal wall. For example, firstly a part of the bladder contour is selected manually as the zone of interests for computing the shear between the bladder and the vaginal. Then $M$ discrete points $\{P_i^0\}$ are sampled on this zone. Given that the shear deformation is initially null, these points are projected orthogonally onto the contour of vagina to obtain $M$ points $\{Q_i^0\}$ who form segments with $\{P_i^0\}$. Then in the presence of shear motion, the segments that are mostly normal to the vaginal walls will deviate from this direction and illustrate the shear effect. We use the variation of angles between the segments as a metric for quantification of the shear strain between organs. The positions of $P_i^t$ and $Q_i^t$ can be calculated by their displacements $dp_i^t$. Let $\vec{t}_i^t$ be the tangent vector of vagina at the point $Q_i^t$ at instant $t$. The shear strain at this position is given by:

$$ \alpha_i^t = \arccos \left( \frac{\overrightarrow{Q_i^t P_i^*} \cdot \overrightarrow{t_i^t}}{\| \overrightarrow{Q_i^t P_i^*} \| \cdot \| \overrightarrow{t_i^t} \|} \right) \quad (2.27) $$

$$ \text{shear}_i^t = \tan \left( \frac{\pi}{2} - \alpha_i^t \right) $$
Finally, after computing the shear strain for all the positions in the discrete domain, we obtain the set of the scalar values representing the shear strain field at instant t (see Figure 2.28 for illustration). Moreover, by combining the spatial distribution of shear strain over time (the dynamic MRI sequence), we will be able to complete a spatio-temporal representation of the shear between organs which is useful for the analysis of the impact of pathologies related to the weakness in the suspension devices.

**Experimental Validation** In this part, we would like to validate our approach with the material points (features) by the comparison of the two methods (the first one being using a classic B-spline). The idea is to apply two methods on the same data set for a quantitative comparison. The captured images consist of an orange, on which a rigid body displacement in the range of 0 – 90 mm was applied. The image spatial resolution is 0.73 mm/pixel and the size is 149×261 pixels. The orange was rotated numerically by an angle of 0 – 20 degrees in order to illustrate the sliding problem because generally the tracking algorithms are very sensible to the rotation. The computed angles of rotation are shown in Figure 2.29. Normally, the angles should be constant along the contour and equal to the imposed values. Nevertheless, it can be observed that the B-spline method is really rotational invariant. Especially, the rotation is not captured at all on the flat portion of the contour. The Bézier method compute better the rotation despite of few slight drifting (sliding effect) also on the flat portion, due to the limitations of the optical flow. As consequence, the values calculated by the Bézier method are much closer to the imposed ones. In order to highlight the difference, we show a more quantitative evaluation in Figure 2.30. The mean angles of rotation computed by the Bézier method remain proportional to the imposed ones. And their standard deviations are less than 1 degree. On the contrary, the result given by B-spline method deviates a lot from the theoretical result.

For the second step, we also compared these two methods for a real case in the context of normal mobility. The objective is to illustrate the drifting effect and then the uncertainty of the computed displacement and shear strain. The test was run on a dynamic MR sequence of a healthy witness which contains 35 images. The spatial resolution is 1.17 mm/pixel and the image size is 256 × 256 pixels. The visual comparison of the results given by two methods is shown in Figure 2.31. The displacements and shear at the middle of straining (the 20th image in the sequence) are used to highlight the difference and the drifting effect. Concerning the result given by the B-spline method, some displacements extremely towards the end of vagina are not reasonable. From the point of view of the physio-anatomy, there exits
<figure><img src="image_65.png" /><figcaption>Fig. 2.29: Comparison of the angles of rotation computed by two methods: the results given by Bézier method are shown on the top and the ones given by B-spline method are on the bottom. The white contour is the initial one. The imposed rotations are 2.5, 10.0 and 17.5 degrees. The color map illustrate the variation along the contour in two cases.</figcaption></figure>
<figure><img src="image_66.png" /><figcaption>**Fig. 2.30:** Comparison of the angles of rotation computed by two methods: the results given by Bézier method are shown on the top and the ones given by B-spline method are on the bottom. The white contour is the initial one. The imposed rotations are 2.5, 10.0 and 17.5 degrees. The color map illustrate the variation along the contour in two cases.</figcaption></figure>

some structures such as the fascias in the zone between bladder and vagina. The bladder, vagina and fascias adhere to each other. Such material is not supposed to stretch considerably without constraint even during the deformation, especially in the normal mobility context. So that the organ walls should neither move freely. Hence the abnormal discontinuity of the displacements on the anterior vaginal wall seems unphysical. The same observation can be found on the posterior vaginal wall. The quantitative comparison is shown in Figure 2.32 As proved in [Lec+15], the zone around the pubis and rectosacral attachment is the low mobility area. From the point of view of anatomy, the bladder and the rectum are attached to these bones during the deformation over time. The two methods give similar results. However, with the B-spline method, we predefined a fix point for each organ as the beginning and ending point (the same point) of the B-spline to position the initial contours. Hence it is normal for the B-spline method to find these low mobility areas. On the anterior and posterior vaginal wall (zone between the vagina and the rectum) and the lower part of the rectum, high mobility can be observed, where the displacements given by the Bézier method attain a magnitude of around 22 mm. Compared to the results in [Lec+15], the displacements are still a bit larger on the anterior vaginal wall. However, through all these comparisons between the results given by the B-spline method, Bézier method and [Lec+15], we are biased in favor of the Bézier method which is more physical.

We performed the same Bézier registration method on a dynamic MR sequence of a patient with prolapse. Because of large displacements in the context of hyper-mobility, the algorithm had some limitations especially when the uterus descended
<figure><img src="image_67.png" /><figcaption>**Fig. 2.31:** Illustration of visual comparison between two method. On the top: displacements at the middle of straining computed by (A) Bézier method and (B) B-spline method. The white contours are the initial ones. On the bottom: (C) initial position of segments illustrating the absence of shear, and deviation of the segments at the middle of straining computed by (D) Bézier method and (E) B-spline method.</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>Fig. 2.32: Quantitative comparison of the displacements computed by two methods: (A,B) Displacement computed by B-spline method at the middle and end of straining, (C,D) Displacement computed by Bézier method at the middle and end of straining.</figcaption></figure>

to the bottom. We would like also mention other limitations due to unwanted effect in the image: the injected gel can flow through the cavity while the patient is pushing and some intrinsic material of the organs can blur the image. All these effect can make the contours deviate from the truly anatomical contours. The restricted sequence from the beginning to the middle of straining contains 14 images, the spatial resolution of which is 1.15 mm/pixel and the size is 240 × 240 pixels. It can be observed that in the high mobility area (posterior vagina) the displacements (Figure 2.33B) are higher than in the normal case (Figure2.32C). The magnitude of displacement is close to 25 mm compared to 15 mm.

The same process was applied on a data set of patient with endometriosis. The sequence contains 32 images of 256 × 256 pixels and the spatial resolution is 1.17 mm/pixel. Due to the presence of wrinkles of the rectum in the image, only the inferior part of the rectum was taken into account for the analysis. In the high mobility area, the magnitude of displacement attained around 10 mm even at the end of straining, which corresponds to the context of hypo-mobility.

**Conclusion** In this section we present a deformable registration approach based on analytical models for multi-organ motion tracking in dynamic MRI sequence. The measured motion is then used for the evaluation of the mobility of pelvic system. The organ was defined by an analytical and parametric model, a composite Bézier curve joined at features (considered as material points). This definition can ensure better the plausibility of the measured motion. The model was then deformed to fit the organ in the dynamic image sequence by using the optical flow method and the virtual image correlation. Finally thanks to the analytical formulation, we computed
<figure><img src="image_69.png" /><figcaption>**Fig. 2.33:** Quantitative evaluation of the displacement field in the context of hyper-mobility: (A) Displacements at the beginning of straining and (B) at the middle of straining. The white contours are the initial ones.</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>**Fig. 2.34:** Quantitative evaluation of the displacement field in the context of hypo-mobility: (A) Displacements at the beginning of straining, (B) at the middle of straining and (C) at the end of straining. The white contours are the initial ones.</figcaption></figure>
the displacement field and shear strain of the moving organs. Hence, the developed
numerical tool, linked with the B-spline based detection, can provide a quantitative
evaluation of organ motion for patient-specific cases. The first perspective of this
work is to assess and validate the computation of the shear strain deformations using
more ground truth data. The quantification of the shear strain could be used to
analyze the abnormal pelvic mobility and to localize the pathological areas related
to the suspension devices which can not be observed on medical images. Such
evaluation would provide objective diagnoses and help in early diagnosis of pelvic
organ prolapse.
Medical and Robotic
Applications

<figure><img src="image_71.png" /></figure>

[...]but I tell you that all these things—yes, from
that star that has just shone out in the sky to the
solid ground beneath our feet—I say that all
these are but dreams and shadows; the shadows
that hide the real world from our eyes. There is a
real world, but it is beyond this glamour and this
vision, beyond these 'chases in Arras, dreams in a
career,' beyond them all as beyond a veil.

— Arthur Machen
The Great God Pan

3.1 Augmented Reality for Per-Operative Guidance

As exposed in the introduction, minimally invasive procedures are quite complex for a surgical standpoint mainly because the visual feedback is relatively incomplete or poor and due to the loss of direct organ manipulation. With recent Augmented Reality techniques, those issues can be circumvented. Internal structures computed from pre-operative scans such as tumors and vessels can be superimposed onto the intra-operative images in order to guide the surgeons during the procedure. Our study is focused on minimally invasive hepatic surgery. In order to establish a full augmented reality system for MIS a number of difficult problems have to be solved [Nic+11]. The abdominal cavity undergoing laparoscopic surgery remains a very challenging environment for computer vision tasks. The surgical instruments interacting with the liver may cause large occlusions, the illumination variations caused by the endoscopic light and the liver bleeding or smokes due to electrocautery may disturb organ motion tracking and 3D structures recovery. Among the most difficult issues is the capacity to tackle, in real time, elastic deformations of the liver. The liver can undergo large deformations due to surgical tools interaction or due to respiratory motion and heart beating. Currently, most existing AR systems
handle rigid motions of organs and only a very limited number of papers address
the problem of elastic organs deformation.

Though the problem of elastic registration or tracking has been largely addressed
in the medical community, elastic deformation in AR systems is only addressed,
to our knowledge, in a very limited number of papers, and mainly applied to 3
organs: kidney, heart and liver. Kidneys do not actually undergo elastic deformations
during surgery but rather move independently from their surrounding tissues. As
such, related works resort to performing a rigid 3D-3D registration e.g. based on
implanted markers [Teb+09] or an Iterative Closest Point (ICP) surface match-
ing [Su+09]. Recently Puerto-Souza et al. [PM13] proposed a matching algorithm
called Hierarchical Multi-Affine capable of long-term tracking for augmentation
during partial nephrectomy. This approach estimates a set of affine transformations
from clusters of features in order to relocate occluded or missing features for a
coherent spatio-temporal mesh registration. Figl et al. [Fig+10] also express the
non-rigid registration problem as a rigid registration between stereoscopic views
and a 4D model of the heart built from 4D scan data, thus leveraging the cyclicity
of the organ motion. A 4D scan of the heart is also used in [Pra+10], but is cou-
pled with a biomechanical model to relate the surface motion to external forces. A
local tuning of the deformation is therefore possible and the surface deformation is
correctly propagated to in-depth invisible structures. No periodic or rigid pattern
can be relied upon in the liver case. Unpredictable, possibly large, deformations can
be observed. Furthermore, volume tracking is required since we are interested in
in-depth augmentation whereas many approaches are dedicated to the tracking of
surfaces. In a similar context, Kim et al. [Kim+12] propose a solution that is both
applicable in a monoscopic setup and integrate a robust image feature reconstruction
approach. But their method only retrieves surface deformation and cannot ensure
an accurate augmentation of inner structures. Mechanical-based deformable models
have proven to be relevant for volume deformation since they allow to define elastic
properties of the shape and thus to infer in-depth structure motion [Spe+11; SLH10;
Ser+03]. Registration is done either by solving mechanical equations given image
constraints provided by reconstructed 3D points [Spe+11] or using the concept of
active models where an energy is minimized which takes into account the internal
behavior of the model -through the biomechanical constraints- and external image
constraints that measure the adequacy between the model and the recovered 3D
data [SLH10; Ser+03].

In this section, we present some results on Augmented Reality for Hepatic surgery
accounting for large liver deformations. We present an approach that relies on the
<figure><img src="image_72.png" /><figcaption>**Fig. 3.1:** Computational flow of the method: The biomechanical model guided by the 3D image-points recovered from intra-abdominal image pair permits to propagate partial tissue deformations to vessels and tumors.</figcaption></figure>

estimation of liver motion to drive a heterogeneous mechanical model capable of deforming the liver. The overview of the method is illustrated in figure 3.1.

### 3.1.1 Mechanical Modeling

#### Parenchyma Model

Most biomechanical studies concerning the constitutive models of the liver parenchyma (see [Ker06] for instance) report the non-linear and viscoelastic behavior of the organ tissue. Since in our case, we do not focus on the transient part of the deformation but rather the static equilibrium under some specific loading conditions, we do not take into account the viscous properties of the tissue. On the other hand, we aim
at modeling large deformations correctly, since during the surgical interventions, important displacements of tissue (e.g. the liver lobes) occur due to the action of surgical tools. For this reason we have opted for finite element method based on co-rotational formulation presented in [Fel00] which allows for large displacements while relying on a linear expression of the stress-strain relationship. Different methods exist for computing the local rotation of each element; in this paper we use a geometric approach proposed in [NPF05] which respects the rotational invariance needed for correct rendering of large deformations (it also has been used in 1.5.1)

In the co-rotational approach, the stiffness matrix K depends on the deformation u and the equation relating the external forces to the displacements can be written as

$$f = K(u)u \quad \text{with} \quad u = x - x_0 \qquad (3.1)$$

where $x_0$ and $x$ represent nodal positions in rest and actual positions, respectively, and $f$ are the external forces.

Assuming that linear tetrahedral P1 elements are employed in the FE formulation of the parenchyma model and the mesh is composed of $N_P$ nodes, the resulting system has $3N_P$ degrees of freedom, i. e. $u_P$, $f_P$ are vectors of size $3N_P$ whereas $K_P$ is a $3N_P \times 3N_P$ matrix, where the subscript $P$ denotes the parenchyma.

It should be noted that our method is not limited to the co-rotational model, which can be replaced by hyperelastic formulation such as the technique based on *multi-plicative Jacobian energy decomposition* presented in [Mar+10].

## Vessel model

The vascular system is regarded as the main source of heterogeneity which has a global influence on the mechanical response of the vascularized tissue. The model employed here is based on work presented in [PDC12]. Besides describing the model in the actual scenario, additional details concerning the assembling of the composite system are given in the following text. Is should be emphasized that a potential viscoelastic response due to the liquid (i. e. blood) circulating in the vessels is not taken into account as only quasi-static scenario is assumed.

From the mechanical point of view, the vascular system is modeled with serially linked beam elements in a similar way as proposed previously in the manuscript (see 1.3) for simulating catheters and guide-wires. We introduce some modifications to the model to take into account the particular nature of vessels, in particular through specific cross section profiles and moments of inertia (see [Prz68] for
details). The static formulation for the deformation of a beam is described by a
system similar to Eq. 3.1 with constituents $u_V$, $f_V$ and $K_V$. However, as each node
is equipped with 6 degrees of freedom due to the rotational components, the size
of vectors $u_V$ and $f_V$ is $6N_V$ and $K_V$ is a $6N_V \times 6N_V$ matrix where $N_V$ is the size of
beam mesh representing the geometry of vessels.

Mechanical coupling

**Mapping of tumors and parenchyma** The size of tumors being relatively small comparing with the whole parenchyma, we assume that its influence on the overall mechanical behaviour is negligible and therefore the coupling with parenchyma can be only geometric. However, a purely geometric mapping may deform the shape of the tumor which may cause erroneous measurements. For this reason the tumor is modeled as a solid object with a very stiff Young's Modulus. This modeling does not affect the performance of the simulation and permits a conservation of the shape of the tumor while allowing correct in-depth propagation of the deformation. With the same formulation as the parenchyma (co-rotational model), coupling between tumors and parenchyma is done using mapping mechanism provided by SOFA [Fau+12] which corresponds to the principle of virtual works used in solid mechanics.

**Coupling between vessels and parenchyma** In order to build the composite model of vascularized tissue, we adopt the method proposed in [PDC12]. Since no relative motion between the vessels and parenchyma is observed in reality, the mapping between the two components can be modeled as a constant coupling. In each step of the simulation, the actual displacements of the parenchyma nodes are mapped to the vessel nodes and reciprocally, the force contribution due to the deformation of the vessel is propagated to the parenchyma. The mapping of forces is based on a principle of virtual work. Additionally, the mechanical contribution of the vessels is added to the global system.

The complete mechanical system of the liver is then: a 3D co-rotational FEM mesh for the parenchyma and the tumors and beam model for the vasculature inside the parenchyma. These 3 models are coupled mechanically (see figure 3.2).
<figure><img src="image_73.png" /><figcaption>Fig. 3.2: 3D heterogeneous biomechanical model of the liver with : (a) heterogeneous liver including the vascular network in wire-frame, in (b) beams generated along the vessels.</figcaption></figure>

### 3.1.2 Non-rigid Registration

Non-rigid registration consists of a 3D point-to-point registration that links the image-based 3D features representing the liver motion to the degrees of freedom of the liver physical model. Registration is done at each frame captured by the endoscopic camera in order to correctly pilot the mechanical model, and thus to estimate the correct position of the underlying tumors. An important key for an efficient registration is the initialization step as it will impact the estimated deformation of the parenchyma and the estimated location of tumors.

#### Initialization

Initialization refers to the step where the alignment of the 3D surface mesh on the laparoscopic image at $t = 0$. This step is critical and is known as an ill-posed problem where we aim at registering two sets of 3D points without knowing the matching between these two sets. In addition, since gas is insufflated (surgical step named as *pneumoperitoneum*) to increase the working space [Ban+12], the preoperative data may no longer correspond to the intra-operative image. Due to this complexity and despite the numerous AR techniques in surgery, very few works have investigated the initial alignment between laparoscopic images and the three-dimensional model. Some have investigated this issue but mostly relies on intra-operative data that are complex to acquire in a clinical routine: for instance [DBB08] use an ultrasound probe to ease the automatic registration process while [Okt+13] use intra-operative
<figure><img src="image_74.png" /><figcaption>**Fig. 3.3:** Initial model alignment (simulated data): The mesh in purple represent the 3D model computed from CT scans in a different reference, the mesh in blue represents the MLS surface reconstructed from laparoscopic image pair. It represents around 30% of the whole liver that is represented in green.</figcaption></figure>

CT scans. Currently, no technique is mature enough to permit an automatic initial alignment with conditions acceptable to clinical routine, and this issue remains a challenge to be solved by the scientific community. This is why, we rely on a manual initial alignment to reduce errors that may emanate from a bad initialization (see figure 3.3): 3D surface reconstruction is built on the first pair of laparoscopic images using the MLS surface approximation and the liver mesh is rigidly aligned using salient geometrical landmarks such as liver contours or surrounded ligaments which enables to provide the rigid transformation (rotation and translation) to align the two datasets. However, this rigid transformation is insufficient to correctly initialize the registration. Indeed, due to gas insufflation, pre-operative data (liver mesh computed for CT-scans) geometrically differs form intra-operative data. To solve this issue, we constrain the biomechanical model (including the vessels and the tumors) with external forces so that it better fits the visible liver surface recovered.

Runtime point correspondence

**Tracking the liver tissue** Tracking object of interest is the most important compo-
nent within Augmented Reality applications. Several previous works have inves-
tigated the motion estimation of non-rigid objects in 3D such as [SF10; WLS12;
Mor+09] but are not applicable in our clinical context because they either rely on
additional data or make strong assumptions on the range of possible deformations.

The visual tracking we used follows the work by Elhawary and Popovic [EP11] who relies on the combination of Lucas-Kanade (LK) [LK+81] optical flow and SURF detector [Bay+08] for robotic-guided endoscopy. This study showed the robustness and accuracy of this combination and its usability for conventional laparoscopic surgery (see figure 3.4 for some experiments). However one major drawback of
<figure><img src="image_75.png" /><figcaption>**Fig. 3.4:** A SURF feature being tracked over time (600 frames) using LK optical flow. (top) tracking performed in both a forward and a backward time direction. Graph shows X and Y coordinates as a function of time. (bottom) 50 × 50 pixel window centered around the feature location at every 100 frames.</figcaption></figure>

this combination is the fact that features loss can happen with the optical flow, preventing to perform long-term tracking. In order to ensure a long term stability for the non-rigid registration, we add an additional layer to the direct tracking. This layer called Cluster-based Filter (CbF) allows to create a set of *n* 3D control points, denoted by a 3Ön coordinate vector *y* that guides the biomechanical model. Instead of constraining the model directly with the extracted features, the CbF computes a displacement field for each control point based on the adjacent features (see figure 3.5). Its benefits are twofold: a) It ensures the minimization to remain stable since the number of control points (and the amount of external forces) is less subject to lost features. Moreover, it permits to avoid to over-constrain the corresponding tetrahedral element which can often tend towards instability. b) it allows to keep only the needed features. Indeed, we can here exploit the capability of the mechanical model to estimate deformations with a reduced number of external loads. The figure 3.6 illustrates the complete approach were a silicone cube is deformed.

**Point correspondence** The temporal non-rigid registration is a point-to-point reg-
istration where we need to map 2D image pixels onto the 3D surface mesh. This
procedure is done using ray casting [Rot82] which is simple to implement and
because it provides several benefits: a) it permits to directly correspond the set of
points to the degrees of freedom of the mechanical model, b) only the 3D points
<figure><img src="image_76.png" /><figcaption>**Fig. 3.5:** Clustering phase: (left) The two views of 3D point cloud estimated from stereo matching. The Blue sphere represents the rejected features. (right) Image is a representation of the clusters. Red sphere represents the Control Points and the Green sphere the essential features (neighbors).</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>**Fig. 3.6:** A sequence of images showing a cube being deformed. (Top) augmented reality images where the mesh in wire-frame is superimposed on the video stream, (Middle) the deformed mesh, (Bottom) the volumetric mesh composed of tetrahedra.</figcaption></figure>
<figure><img src="image_78.png" /><figcaption>Fig. 3.7: Initial model alignment (simulated data): The mesh in purple represent the 3D model computed from CT scans in a different reference, the mesh in blue represents the MLS surface reconstructed from laparoscopic image pair. It represents around 30% of the whole liver that is represented in green.</figcaption></figure>

that intersect the liver surface after ray-casting are kept, the features that do not belong to the liver are filtered out from laparoscopic images (see figure 3.7).

### 3.1.3 Implementation details

Non-rigid registration can be seen as a stretching energy minimization between the three dimensional features recovered from laparoscopic images that represent the tracking energy and the biomechanical model derived from preoperative CT data. External forces are defined by pairing the m 3D points $y$ to the n degrees of freedom $x$ of the biomechanical model. $f_t(y) = k(y - y^0)$ where $y^0$ defines the tracked control points at time $t = 0$ and $k$ can be defined as a stiffness. The tracked coordinates can be expressed using barycentric coordinates of adjacent degrees of freedom $y = Lx$.

To be as accurate as possible, boundary conditions representing the falciform ligament (ligament that attaches the liver to the front body wall) and the underlying fat needs to be defined. These boundaries conditions represents fixed parts of the liver that could not move or deform over time. The figure 3.8 illustrates boundary conditions on our model.
<figure><img src="image_79.png" /><figcaption>Fig. 3.8: Boundary conditions: The liver model (in green) is constrained by fixed vertices (in purple) that represent the falciform ligament and the underlying fat.</figcaption></figure>

Given these additional concerns, the simulation is then performed at every time steps and external forces due to image tracking as well as internal forces due to liver model are computed.

### 3.1.4 Experimental Results

The validation of such framework is complex as clinical context makes it difficult to build an assessable, reproducible set-up in order to evaluate the performance of each element of the proposed pipeline. The thesis of Nazim Haouchine [Hao15] details the methodology we used in order to validate our approach: virtual validation with computer-generated data, in-silico experiment with silicone phantom, ex-vivo experiments with pig liver and in-vivo experiments with human data. The figure 3.9 illustrates some images extracted from a AR sequence where the liver has been successfully deformed based on the motion tracking. For the complete study and the quantitative results, we encourage the reader to refer to the PhD thesis of Nazim Haouchine.

## 3.2 Inverse FEM simulation for adaptive radiotherapy

Radiation therapy (or *radiotherapy*) is one of the possible treatments for head and neck cancers. It uses high-energy X-rays to destroy the cancer cells. A treatment is
<figure><img src="image_80.png" /><figcaption>**Fig. 3.9:** A sequence of images showing the superimposition of the real-time biomechanical model onto the human liver undergoing deformation due to surgical instrument interaction. On the left virtual camera with the original laparoscopic angle and right a different angle.</figcaption></figure>
established by using a treatment planning system (TPS) [Web03] which combines patient medical images, radiation transport simulations and optimization algorithms in order to expose tumors to X-rays while sparing healthy structures. The treatment plan is then applied 5 days per week during 6 to 7 weeks in order to destroy the tumors. During these 7 weeks, the patient is exposed to several side effects, and in particular an important weight loss. This induces the motion and deformation on the anatomical structures and the TPS does not take into account these topography changes, which may lead to an important X-ray exposure of healthy tissues [Nel+12] [Vei+14]. For instance it is reported [DNV06] that xerostomia (loss of saliva) is due to a higher (than planned) exposure of the parotid glands while treating throat cancers.

The main objective in this section is to adapt the planning to account for morphological modifications in order to limit the radiation exposure of healthy structures. It has been shown that non-rigid registration and daily computation of the dose can reduce the radiations [Vei+14]. But the challenge remain on the registration method over the 7-week period. While significant works have been achieved recently in the field of automatic non-rigid registration (the reader may refer to [CHH04] for a recent survey), these methods do not provide an easy control for the physicians. These algorithms also lack robustness and consistency when images are very complex. On the contrary, dealing with manual segmentations and/or registrations is time-consuming for the physicians and is not a viable solution for adapting the planning along with the treatment of the patients.

Our approach will be the same as in the previous section: using images to guide FEM models that represent the anatomical structures we want to track over time. But we want to improve it by avoiding to apply forces (or displacements) through the analysis of images since the physics of these boundary conditions is not correct (for a mechanical standpoint) as it mixes physics-based forces with image-based forces. As a result we will solve an inverse FEM simulation where several landmarks will guide the resolution of the problem.

### 3.2.1 Real-time Inverse FEM

Starting from what has been described in section 1.2, we extend the equations used for the resolution of real-time static FEM model. During each step i of the simulation, a linearization of the internal forces is computed:

$$f(x_i) \approx f(x_{i-1}) + K(x_{i-1})dx \quad (3.2)$$
where f provides the volumetric internal stiffness forces at a given position x of the nodes, K(x) is the tangent stiffness matrix that depends on the actual position of the nodes and dx is the difference between consecutive positions in time dx = xᵢ - xᵢ-₁. The lines and columns that correspond to fixed nodes are removed from the system to improve the condition number of the matrix K. Static equilibrium (the sum of external and internal force equals to zero) is sought at each step:

$$-K(x_{i-1})dx = p + f(x_{i-1}) + J^T \lambda \quad (3.3)$$

where p represents the external forces (e.g. gravity) that are known and J^Tλ gathers the contributions of the Lagrange multipliers. Three types of multipliers are defined:

* **boundary multipliers** λ_b: we use these constraints to describe the external efforts applied on the boundary conditions that creates the deformation. The location of the boundary conditions are supposed to be known, the directions of the effort J^T can be updated at each step, and λ_b is the unknown intensity of the efforts on boundaries. We can set (and update at each step i) an interval of potential values min ≤ λ_b ≤ max.

* **parameters multipliers** λ_p: these parameters influence the computation of the internal forces. We use a local derivation of the internal force by the parameter p: f(x, p + dp) ≈ f(x, p) + (δf(x, p)/δp)dp. In that case, J^T = δf(x, p)/δp and λ_b = dp is the variation of the parameter. To keep the validity of the local derivation over a step i, we can set -ε ≤ λ_p ≤ ε.

* **registration multipliers** λ_r: set interactively by the user to do a local manual registration on a small number of points.

Contrary to a lot of existing registration methods, we do not put any force (or energy) to the association of points or on image pixels so λ_r = 0. Even if null, these multipliers are useful to build the optimization problem.

Indeed, the next step consists of the projection of the FEM model equations into the constraint space: the size of matrix K is often very large so an optimization in the motion space would be computationally very expensive. Instead, using the Schur complement (also called Delassus operator) of the constraint problem, we do a projection that dramatically reduces the size of the research space.

Three steps are followed, that are standard in a constraint solving process:
* **Step I**, a free configuration $x_{\text{free}}$ of the deformable model is found by solving equation (3.3) with $\lambda = 0$. For the constraint defined on registration point, we compute a violation noted $\delta_r^{\text{free}}$ which provides a vector between the registered position of the points and the position given during the free motion.

* **Step II:** This step is central in the method. It consists in projecting the mechanics into the constraint space. As the constraints are the inputs (registered points) and outputs (parameters and efforts on boundary) of the inverse problem, we obtain the smallest possible projection space for the inverse problem:

$$ \delta_r = \underbrace{[\mathbf{J}_r \mathbf{K}^{-1} \mathbf{J}_p^T]}_{\mathbf{W}_{rp}} \lambda_p + \underbrace{[\mathbf{J}_r \mathbf{K}^{-1} \mathbf{J}_b^T]}_{\mathbf{W}_{rb}} \lambda_b + \delta_r^{\text{free}} \quad (3.4) $$

$\delta_r$ represents a vector between registered and actual positions of points chosen by the user. Then a Quadratic Programming (QP) problem is set by minimizing the norm of this vector.

$$ \min(\frac{1}{2}\delta_r^T \delta_r) = \min(\frac{1}{2} \begin{bmatrix} \lambda_p \\ \lambda_b \end{bmatrix}^T \begin{bmatrix} \mathbf{W}_{rp}^T \\ \mathbf{W}_{rb}^T \end{bmatrix} [\mathbf{W}_{rp} \quad \mathbf{W}_{rb}] \begin{bmatrix} \lambda_p \\ \lambda_b \end{bmatrix} + \begin{bmatrix} \lambda_p \\ \lambda_b \end{bmatrix}^T \begin{bmatrix} \mathbf{W}_{rp} \\ \mathbf{W}_{rb} \end{bmatrix} \delta_r^{\text{free}}) $$

$$ \text{subject to } min \le \lambda_b \le max \quad \text{and} \quad -\epsilon \le \lambda_p \le \epsilon \quad (3.5) $$

The size of the QP problem is much smaller than solving the problem in the motion space of equation 3.3, allowing to solve this problem in real-time. In practice, the QP matrix is always positive and definite iff the size of $\delta_r$ is greater than the number of optimized values in $\lambda_p$ and $\lambda_b$.

* During **Step III**, the final configuration of the deformable model is corrected by using the value of the constraint response using $x_i = x_{\text{free}} + \mathbf{K}^{-1}(\mathbf{J}_p^T \lambda_p + \mathbf{J}_b^T \lambda_b)$. In practice, we use a $LDL^T$ factorization of the matrix $\mathbf{K}$ and not $\mathbf{K}^{-1}$ during the computation.

### 3.2.2 Experimental validation

In this section, we present the application of the method in the context of radiotherapy of the neck (throat cancer) where patient weight loss induces deformations of anatomical structures that are not taken into account in the treatment. At the beginning of each therapy session, a scan of the patient is taken and a rigid registration of the planning to the actual position of the patient is realized. In order to facilitate this
registration, the patient is wearing a rigid mask. Notably, the deformations of the parotids make them move towards the center of the neck and eventually intersect the main target volume of the therapy (see figure 3.10). Consequently parotids may be irradiated more than initially planned leading to xerostomia (20% of patients).

<figure><img src="image_81.png" /><figcaption>Fig. 3.10: Volume loss of parotids: (Left) segmentations of the parotids at weeks 1 (red) and 6 (blue). It is worth noticing the volume loss of the parotid as well as the motion of the center of mass. These two parameters have been used to characterize the deformation of parotids in [Bar+04]. (Right) Due to weight loss, parotids may intersect the target volume (in yellow).</figcaption></figure>

Yet, with a robust registration performed just before the therapy, the planning could be adjusted. We emphasize that the parotids are poorly visible in the images, making the automatic registration not robust. With our inverse approach, the radiotherapist can use his expertise (knowledge of anatomy and meaningful parameters used in medical studies) to perform the registration and he/she will have a direct control of the parameters used for the registration. For instance [Bar+04], studied the parotids deformation by performing CT scans of the patient three times a week during the entire treatment. From that study, the deformation of the parotids is characterized by two parameters: the volume loss of the parotids and the motion of its center of mass (due to the volume loss of other structures). Parotids deformations observed are large (more than 30 % the size of the structure) and can not be captured with linear models.

The current medical routine does not adapt the treatment since it involves the manual segmentation of the structures -which is time and manpower consuming- and the computation of the new planning. However our method can dramatically reduce the time required to adapt the planning while achieving comparable accuracy to manual segmentation. We tested our approach on a *ground truth* set of 7 patient datasets that contains the 3D images of the CT scan done every week of the therapy (total: 7*6 images). Comparison between manual segmentations of the parotids (performed by the radiologists) and our method is achieved on all available images (42) by computing the DICE coefficient. A single dataset (6 images) has been manually segmented by two radiologists and an average DICE coefficient of 0.7 has been computed to serve as a reference for the quality of our method. On these data, our method can be executed very quickly (completion of the registration is done in a
single minute) with respect to a full manual segmentation making it compatible with the time constraints of a clinical routine. The graph 3.11 (left) illustrates that the parotids deformation is significant and second that our method exhibits really good similarity compared to manual segmentation (average DICE between [0.8;0.9]).

<figure><img src="image_82.png" /><figcaption>Fig. 3.11: Validation: (left) similarity between the initial segmentation and *ground truth* geometry in blue curve illustrates the deformation of the parotids (DICE decreasing), the red curve exhibits the good similarities between our semi-automatic registration and the *ground truth* geometry; (middle) planning adaptation using our registration vs no planning adaptation (right). The measured radiation is much lower when the planning is adapted.</figcaption></figure>

**Efficiency of the whole approach:** A dataset was selected for which the deformations were important and the parotids were not infiltrated by the tumor (therefore out of the target volume). We have a closer look at the last session of the therapy and particularly at the irradiation map of the parotids without considering planning adaptation 3.11 (right) and with planning adaptation using our method to register the right parotid (middle). The resulting maps from the TPS show that the irradiation of the right parotid is significantly reduced and may limit the appearance of irradiation side-effects.

## 3.3 Motion Control for Interventional Radiology Procedures

Using FEM in an inverse problem as in the previous section has been then an approach that we have used in many robotics application for instance to build a closed-loop control for a robot [Zha+16] or to perform visual servoing [Zha+17] or to combine computer vision and inverse model to use the robot as a cheap force sensor [ZDD18b]. We choose not to detail all these works as the reader can refer to them, but we will focus on a work that combine robotics and interventional radiology. We already have explained in the previous chapters (more specifically in
chapter 1.5.2) that interventional radiology procedures requires excellent eye-hand coordination and strong expertise to anticipate how the radiology devices (catheters, guide-wires...) will behave inside the vascular network. Complementary to the training and planning tools we have proposed, a clinical application providing the automatic or semi-automatic control of a robotized catheter actuated by cables. Without removing the radiologist from the loop, we would like to consider a clinical scenario where the radiologist defines a path that the catheter must follow given some per-operative medical images and our control algorithm actuates the robotized catheter to follow this trajectory. This work has been conducted in collaboration with Junghwan Back and Hongbin Liu from Kings College London and the catheter robot used was an upgraded version of [Bac+18]. This catheter robot has a polymer material structure which is MR-compatible, has a diameter of 7mm in order to

<figure><img src="image_83.png" /><figcaption>Fig. 3.12: Catheter robot from Kings College London with a polymer structure and 4 tendons that enables the robot tip to bend.</figcaption></figure>
### 3.3.1 Simulation with actuation and contacts

The control of soft robot like the catheter we use requires the simulation of soft structures. The simulation of soft robots is no different from the simulation of anatomical structures and as a consequence we will use the methods we have detailed so far in this manuscript. More specifically the control of the robot requires to solve an inverse problem and we will use the methodology described in the previous section. Given the equilibrium equation, we use the Schur complement of the constraint problem, then the FE model equation can be projected into the constraint space that drastically reduces its size. Then, we can compute the following variables at any step k given variables at step k − 1:

$$ \delta_{e,k} = \mathbf{W}_{ea}(\mathbf{x}_{k-1}) \lambda_{a,k} + \mathbf{W}_{ec}(\mathbf{x}_{k-1}) \lambda_{c,k} + \delta_{e,k}^{free} \quad (3.6) $$

$$ \delta_{a,k} = \mathbf{W}_{aa}(\mathbf{x}_{k-1}) \lambda_{a,k} + \mathbf{W}_{ac}(\mathbf{x}_{k-1}) \lambda_{c,k} + \delta_{a,k}^{free} \quad (3.7) $$

$$ \delta_{c,k} = \mathbf{W}_{ca}(\mathbf{x}_{k-1}) \lambda_{a,k} + \mathbf{W}_{cc}(\mathbf{x}_{k-1}) \lambda_{c,k} + \delta_{c,k}^{free} \quad (3.8) $$

where the following variables are defined as:

* subscript a refers to the actuators, subscript e refers to the end-effector and subscript c refers to the contacts

* $\delta_{e,k}$ Position vector of effectors at step k

* $\delta_{a,k}$ Position vector in the actuator space at step k

* $\delta_{c,k}$ Position vector in the contact space at step k

* each matrix $\mathbf{W}_{ij}$ is defined as $\mathbf{W}_{ij}(\mathbf{x}_k) = \mathbf{J}_i(\mathbf{x}_k) \mathbf{K}^{-1}(\mathbf{x}_k) \mathbf{J}_j^T(\mathbf{x}_k)$. The couples $(i, j)$ could take the values $(i, j) = (e, a); (e, c); (a, a); (a, c); (c, a); (c, c)$. Then:

- $\mathbf{W}_{ea}(\mathbf{x})$ is the compliance matrix between effectors and actuators,

- $\mathbf{W}_{ec}(\mathbf{x})$ is the compliance matrix between effectors and contacts,

- $\mathbf{W}_{aa}(\mathbf{x})$ is the compliance matrix between actuators

- $\mathbf{W}_{ac}(\mathbf{x})$ is the compliance matrix between actuators and contacts,
- $W_{ca}(\mathbf{x})$ is the compliance matrix between contacts and actuators,

- $W_{cc}(\mathbf{x})$ is the compliance matrix between contacts.

* $\lambda_{a,k}$ is the force contribution vector of actuators (at step $k$)

* $\lambda_{c,k}$ is the force contribution vector of contacts (at step $k$)

* $\delta_{a,k}$ is the position vector in the actuator space (at step $k$). $\delta_{a,k}^{free}$ represents the corresponding position vector when no forces are applied ($\lambda_{a,k} = 0$ and $\lambda_{c,k} = 0$)

* $\delta_{e,k}$ is the position vector of effectors (at step $k$). $\delta_{e,k}^{free}$ represents the corresponding position vector when no forces are applied ($\lambda_{a,k} = 0$ and $\lambda_{c,k} = 0$)

* $\delta_{c,k}$ is the position vector in the contact space (at step $k$). $\delta_{c,k}^{free}$ represents the corresponding position vector when no forces are applied ($\lambda_{a,k} = 0$ and $\lambda_{c,k} = 0$)

Without contacts and with actuation only, the position of the soft robot is defined as:

$$ \mathbf{x}_{k+1} = \mathbf{K}_k^{-1} \mathbf{H}_a^T \cdot \lambda_k + \mathbf{x}_k^{free} \quad (3.9) $$

When contacts are considered, the collision response is based on the Signorini's law [Sig33]. For every contact point, there is a complementarity relation between the interpenetration gap $\delta_{c,k}$ and the contact force $\lambda_{c,k}$ among the normal direction. The complementarity relation can be written as [DAK04]:

$$ 0 \le \delta_{c,k} \perp \lambda_{c,k} \ge 0 \quad (3.10) $$

where $\perp$ is the complementarity operator. It states that each step $k$ of the simulation either the forces due to contacts ($\lambda_{c,k}$) are null or that no DOF is in contact with obstacles ($\delta_{c,k} \neq 0$). Using equations 3.8 and 3.10, the LCP could be written as:

$$ \mathbf{W}_{ca}\lambda_a + \mathbf{W}_{cc}\lambda_c + \delta_c^{free} \ge 0 \quad (3.11) $$

$$ \lambda_c \ge 0 \quad (3.12) $$
$$ \lambda_c^T (\mathbf{W}_{ca}\lambda_a + \mathbf{W}_{cc}\lambda_c + \delta_c^{free}) = 0 \qquad (3.13) $$

### 3.3.2 Control design

In a constrained environment, [YC14] have highlighted two canonical problems for the control design of continuum robots if contacts are not taken into account in the closed loop design: artificial singularity and inverted mapping (see figure 3.13). Without taking contacts into accounts, the control could lead to a positive-feedback loop. Both the artificial singularity and the inverted mapping are generated by the coupling between the insertion and bending. Therefore, to avoid the two canonical problems, we propose a decoupled control strategy by decoupling the insertion and the bending. Besides, it aims to increase the insertion safety because it enables automatic bending control and manual insertion control.

<figure><img src="image_84.png" /><figcaption>Fig. 3.13: Two canonical cases described where the contacts highly influence the kinematics of the catheter robot. Arrows are the motion directions of the tip for imposed motion at the base: insertion, right cable movement, left cable movement. (a) artificial singularity and (b) inverted mapping result in the deterioration and even the positive-feedback loop.</figcaption></figure>

The proposed closed-loop strategy (see figure 3.14) needs to combine both the catheter robot and its FE model which is simulated in real-time. The FE model is employed to compute the Jacobian matrix and to detect collisions between the robot
and the environment. The same insertion input $\delta_{in}$ is sent to both the catheter robot and its FE model.

<figure><img src="image_85.png" /><figcaption>Fig. 3.14: Framework of decoupled closed-loop bending control for catheter robot through contacts. For the implementation, the FE model of catheter robot and its environment is simulated in real-time. The insertion and bending are decoupled by the projection block which is introduced in Section IV.A. Two independent controllers are employed for the control of catheter robot and its FE model respectively. $\delta_d$ and $\delta_{be}$ are respectively the desired position and the bending control inputs for the catheter robot, while and $\delta_{FEM}$ are the corresponding variables for the FE model. $\delta_{in}$ insertion input for both the catheter robot and the FE model.</figcaption></figure>

## Projection step

The projection component is used to decouple the insertion and the bending for the catheter robot. The movement of catheter robot is modeled in 3D space defined by the coordinate system OXYZ. To achieve the projection, the local coordinate system Pxyz is computed at each step based on the 3D tip position and the predefined trajectory so that it moves with the insertion of catheter robot. After the projection process, the movement of catheter tip in 3D space is converted to the insertion along the tangent axis x and the bending movement in 2D plane defined by the normal axis z and the projection axis y.

## FEM Controller design

The objective of the closed-loop controller for the FE model is that the effector of the FE model ($\delta_e$) converges towards the desired location ($\delta_d$). We can define an optimization problem where we want to find the forces $\lambda$ that minimizes the objective function $\Gamma$:
$$ \Gamma(\lambda) = \frac{1}{2} \| \mathbf{W}_e \lambda + \delta_e^{\text{free}} - \delta_d \|^2 \quad (3.14) $$

We can add some constraints on the forces that should be generated such as $0 \le \lambda_c \le \lambda_{c,max}$ to limit the magnitude of contact forces or $\lambda_{a,min} \le \lambda_a \le \lambda_{a,max}$ to take into account the physical limitations of actuators. This optimization problem is rewritten as a QP-problem that can be solved using dedicated software packages. The reader may refer to our journal paper [Zha+19b] for equations details. In this end we obtain a standard QP-problem:

$$ \begin{array}{ll} \min_{\lambda} & \frac{1}{2} \lambda^T \mathbf{Q} \lambda + \mathbf{c}^T \lambda \\ \text{s.t.} & 0 \le \lambda \le \lambda_{\max} \\ \text{and} & \mathbf{W}_c \lambda + \delta_c^{\text{free}} \ge 0 \end{array} \quad (3.15) $$

where we make the assumption that the defined trajectory is always reachable. In practice, with assumption of static environment, we plan it off-line and make sure that it is always reachable by our robot.

## Robot Controller design

Given equations formulated in subsection 3.3.1, we can express the position of the effector given compliance matrices in:

$$ \delta_{e,k+1} = \delta_{e,k} + \mathbf{J}(\mathbf{x}_k) \Delta \delta_{a,k} \quad (3.16) $$

where $\mathbf{J}(\mathbf{x}_k) = \mathbf{W}_{ea}(\mathbf{x}_k) \mathbf{W}_{aa}^{-1}(\mathbf{x}_k)$ is the Jacobian matrix of the robot. Given this Jacobian, the control vector $U_k = \mathbf{J}(\mathbf{x}_k) \Delta \delta_{a,k}$ allows to rewrite the kinematic model as:

$$ \delta_{e,k+1} = \delta_{e,k} + U_k \quad (3.17) $$

The tracking error is defined as $e_k = \delta_{e,k} - \delta_{d,k}$. In the task space, the control vector $U_k$ can be designed as:

$$ U_k = -K_P e_k - K_D (e_k - e_{k-1}) \quad (3.18) $$
where $K_P$ and $K_D$ are constant parameters for the proportional gain and differential gain respectively. Then, the Pseudo-inverse based control allocation is employed to obtain a unique solution:

$$u = \mathbf{J}^{+}(\mathbf{x}_k) U_k \qquad (3.19)$$

$$u = -\mathbf{J}^{+}(\mathbf{x}_k) [K_P e_k + K_D (e_k - e_{k-1})] \qquad (3.20)$$

### 3.3.3 Experimental validation

In order to validate this approach, we use a commercial electromagnetic tracking system (Aurora tracking from NDI company) with the sampling frequency 20Hz and the sub-millimetric accuracy. A rigid tube (diameter 28mm) is fixed with respect to the base of the catheter robot. Positions of seven feature points are captured by the tracking system and are employed for the calibration between the coordinate system of the sensor and that of the robot. The movement of catheter robot consists of bending and insertion.

<figure><img src="image_86.png" /><figcaption>Fig. 3.15: Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</figcaption></figure>

The positioning results are shown in figure 3.16 which includes two cases with different catheter length (60 mm and 166 mm). The tracking error converges faster when the length of catheter is shorter. For both cases, the positioning errors converge to a small bound (around 0.2mm). Usually a longer length of catheter robot generates a larger convergence bound if there is no interaction with the
environment. However, in our experiments, the contacts limit the deflection and therefore reduce the convergence bound.

<figure><img src="image_87.png" /><figcaption>**Fig. 3.16:** Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</figcaption></figure>

Additional experiments were conducted using only simulation in order to focus on the control algorithms themselves and can get rid of the sensor noise, imprecise modeling of robot and contacts. As shown in figure 3.17, the catheter robot can be controlled from the initial configuration (a) to the final configuration (b) with the decoupled control strategy. If the insertion and bending are coupled between each other, the Jacobian matrix is inverted during the movement. However, the decoupled control strategy can successfully avoid the problem of the Jacobian matrix rotation when the robot interacts with the contacts.

The trajectory tracking can also be achieved when the catheter robot interacts with soft contacts. The decoupled control strategy is employed and validated using simulation which is shown in figure 3.18. To estimate the Jacobian matrix and contact forces, both catheter robot and soft tube are modeled using FE method and are simulated in real-time.

## 3.4 Perspectives

The three different applications presented in this chapter illustrate how the foundational works presented in the manuscript could open up for new high-value applications in different fields. The combination of interactive (or at least *fast*) simulation of
<figure><img src="image_88.png" /><figcaption>**Fig. 3.17:** Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</figcaption></figure>

<figure><img src="image_89.png" /><figcaption>**Fig. 3.18:** Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</figcaption></figure>

deformable objects with computer vision and automatic control could provide many beneficial outcomes in surgery, interventional radiology or robotics. We also want to emphasize on the fact that these applications require real multi-disciplinary approach to be really efficient: some components are developed throughout the manuscript but some technical others like software programming, algorithm optimization and modular software platform are mandatory to achieve these applications. A standard unified platform that includes all these aspects would help the simulation and robotics community but the standard is still to emerge. SOFA<fnref n="1" /> is an attempt to answer these technical deadlocks and all the applications or simulations presented in this manuscript were built using SOFA.

<footnote n="1">¹www.sofa-framework.org</footnote>
Conclusion

” I didn't think; I experimented.

— Anthony Burgess

This manuscript presented several research works conducted these past year. Focussing on interactive simulations for deformable objects, we have proposed some contributions to the modeling of deformable solids, to the optimization of simulation and some training applications that can be derived from these works. We also shown that these simulations need to rely on accurate data to be really useful. Since medical images are used extensively in clinical routine, we have proposed some methods to build simulation-compatible data from images and we have also illustrate how the data, if acquired in run-time, could be used to guide or increase the prediction capabilities of simulation. The combination of interactive simulations and computer vision algorithm have led to several applications whose outcomes are more than training tools... and we have illustrate some examples such as a planning tool for radiotherapy or an augmented reality tool that enables interactive per-operative guidance for surgeons. These successful combination is not limited to augmented reality / virtual reality applications for medicine but can also be employed in robotics. The algorithmic principles can be derived and using automatic control methods we have developed some applications for automatic controls of soft robots.

Several perspectives for my research works can be issued from this summary and the state-of-the-art. A lot of attention has been directed to artificial intelligence / supervised (deep) learning and these past few years have seen numerous papers using model-free methods based on AI to solve simulations or automatic controls issues... this has led to an opposition between model-free researchers and model-base researchers. While the benefits of supervised learning are indisputable, the main limitation is the huge volume of data to train the neural networks used for the learning. Some application contexts can provide relatively easy all the data needed but our experience in medical applications shows that large and most important relevant data are really complex to acquire: medical images are abundantly available but they should pass through several filters some social (ethical committees, anonymization...) some technical (segmentation, reconstruction...) and there are still some limitations such as the acquisition speed in case of dynamic images or the
fact that these images contain mostly geometric data... which are not sufficient for
mechanical simulation. Our vision for the next few years to combine the prediction
capabilities of AI with the scenario generation of simulations: simulation can create
a virtual testbed that can produce / reproduce many clinical procedures, pathologies,
complications... that can be recorded to be used to train the supervised learning
algorithms. For that context, interactive applications is not necessarily a requirement
but the faster the simulations the faster it can produce the amount of data needed.
Thus all the deadlocks related to simulation speed or simulation optimizations that
we have mentioned in chapter *Interactive Physics-based Simulation* are still relevant
and worthy to investigate. We had a collaboration with a company (Robocath) about
the automatic navigation of catheter in cardiac interventional radiology that uses
the combination of approaches and the results are promising. Some of these results
have been published in [Sch+22] and [Zha+22] and we plan to extend them in the
future.

A second perspective is mostly related to the previous one, as (dynamic) 2D / 3D
imaging techniques are available through medical imaging devices, 3D scanners or
cameras... the amount of geometric data generated is huge. These data, while di-
rectly usable for visualization or rendering, are not suited for mechanical simulation:
topologies may exhibit some unexpected artifacts like unwanted holes, very low
mesh quality (lots of degenerated triangles)... As a consequence, a time-consuming
and tedious post-processing step is mandatory and most of the times this is done
with the expertise of someone that possess simulation and geometric modeling skills.
Moreover since modeling softwares uses triangular mesh representation for virtual
objects, they lack the ability to define volumetric properties that are required by 3D
mechanical simulation and they are the best representation for 3D objects. In the
chapter *Geometric Modeling for Interactive Simulation*, we have made an attempt
to use implicit surfaces as object representation because of its intrinsic properties
(geometric continuity, smoothness) and because it speeds-up drastically the compu-
tation time of the collision detection step in the simulation. With the rise of new
graphics cards that include hardware ray casting (which can be used for collision
detection as well as rendering) like in the RDNA Architecture from AMD, we could
expect significant speed-ups in the simulation by using these features with implicit
surfaces. The major scientific deadlock remains how to build implicit surfaces with
an arbitrary geometry (be it a triangular mesh or a point cloud). What we have
proposed in the second chapter is relevant for tubular structures but it cannot be
derived for arbitrary shapes. Therefore we want to investigate this field with having
always interactive simulation as a target (thus aiming for a compact representation
of surfaces).
A third perspective is related to the difficulties encountered for the experimental validation of our approaches. We lack ground truth data in the medical simulation community or in the soft robotics community. Since we aim at providing methods that can be derived into medical applications used in clinical routine or robotics applications used in industry or surgery, special care should be taken to validate experimentally in scenario that are close to *in-vivo* situations. A gradual approach ranging from computer validations to *in-vivo* validations stepping by *in-silico* and *ex-silico* validations appears to be relevant, it cannot be done in every context. Silicone phantoms that we have used numerous time in the chapter *Medical and Robotic Applications* enables to provide a reproducible and assessable testbed for validations. They, however, lack dynamic features (motion and deformation) that can mimic anatomical soft structures. Our proposal for the next years is to build controllable and instrumented soft robots that reproduce anatomical soft bodies. This research direction combines simulation, automatic control and soft robots design in order to produce a validated testbed where motion can be controlled (to mimic cardiac cycle for cycle), elastic deformations are known or measurable, forces are measured... over time. We think this step is mandatory to fasten the time needed to go from the research paper to the pre-industrial prototype. In order to accelerate on this perspective, I submitted a proposal for an european call that includes researchers from university of Copenhagen, Universidad Rey Juan Carlos (Madrid) and Vrije Universiteit Brussel to build soft robotized twins of organs and anatomical structures.

These three research directions are not only extensions of the three chapters pre-
sented in this manuscript. From the capitalization of our research experience these
past 10 years and the integration of recent works and recent collaborations, we
have identified these three topics were we can have an impact and provide original
contributions. The first two can be investigated in parallel while the third one is
more transverse and multi-disciplinary and will require a network of collaborations
to fill the different expertise we need.
<figure><img src="image_90.png" /></figure>
Bibliography

[AC04] Gallagher AG and Cates CU. "Approval of virtual reality training for carotid stenting: What this means for procedural-based medicine". In: *JAMA* 292.24 (2004), pp. 3024–3026 (cit. on p. 2).

[AKN07] T. Alderliesten, M.K. Konings, and W.J. Niessen. "Modeling Friction, Intrinsic Curvature, and Rotation of Guide Wires for Simulation of Minimally Invasive Vascular Interventions". In: *Biomedical Engineering, IEEE Transactions on* 54.1 (Jan. 2007), pp. 29–38 (cit. on p. 54).

[Bac+18] Junghwan Back, Lukas Lindenroth, Kawal Rhode, and Hongbin Liu. "Three dimensional force estimation for steerable catheters through bi-point tracking". In: *Sensors and Actuators A: Physical* 279 (2018), pp. 404–415 (cit. on p. 130).

[Ban+12] Jordan Bano, Alexandre Hostettler, SA Nicolau, Stephane Cotin, Christophe Doignon, HS Wu, MH Huang, Luc Soler, and Jacques Marescaux. "Simulation of pneumoperitoneum for laparoscopic surgery planning". In: *International Conference on Medical Image Computing and Computer-Assisted Intervention*. Springer. 2012, pp. 91–98 (cit. on p. 118).

[Bar+04] Jerry L Barker Jr, Adam S Garden, K Kian Ang, Jennifer C O'Daniel, He Wang, Laurence E Court, William H Morrison, David I Rosenthal, KS Clifford Chao, Susan L Tucker, et al. "Quantification of volumetric and geometric changes occurring during fractionated radiotherapy for head-and-neck cancer using an integrated CT/linear accelerator system". In: *International Journal of Radiation Oncology* *Biology* *Physics* 59.4 (2004), pp. 960–970 (cit. on p. 128).

[Bar+94] Victor Barrios, J Torres, Guillermo Montilla, Lilia Hernandez, Naykiavic Rangel, and A Reigosa. "Cellular edge detection using a trained neural network explorer". In: *Engineering in Medicine and Biology Society, 1994. Engineering Advances: New Opportunities for Biomedical Engineers*. Proceedings of the 16th Annual International Conference of the IEEE. IEEE. 1994, pp. 1075–1076 (cit. on p. 54).

[Bay+08] Herbert Bay, Andreas Ess, Tinne Tuytelaars, and Luc Van Gool. "Speeded-up robust features (SURF)". In: *Computer vision and image understanding* 110.3 (2008), pp. 346–359 (cit. on p. 119).

[BM94] Dennis G Begos and Irvin M Modlin. "Laparoscopic cholecystectomy: from gimmick to gold standard." In: *Journal of clinical gastroenterology* 19.4 (1994), pp. 325–330 (cit. on p. 23).

[Ber+13] Matthew Berger, Joshua A Levine, Luis Gustavo Nonato, Gabriel Taubin, and Claudio T Silva. "A benchmark for surface reconstruction". In: *ACM Transactions on Graphics (TOG)* 32.2 (2013), pp. 1–17 (cit. on p. 63).
[BHC93] JC Bezdek, LO Hall, and LP Clarke. "Review of MR image segmentation techniques using pattern recognition". In: *Medical Physics* 20.4 (1993), p. 4 (cit. on p. 54).

[Bil14] Alexandre Bilger. "Patient-specific biomechanical simulation for deep brain stimulation". PhD thesis. Université des Sciences et Technologie de Lille, 2014.

[BTC95] Eric Bittar, Nicolas Tsingos, and Marie-Paule Cani. "Automatic Reconstruction of Unstructured 3D data: Combining Medial Axis and Implicit Surfaces". Anglais. In: *Eurographics*. Ed. by Frits Post and Martin Göbel. Vol. 14. 3. Published under the name Marie-Paule Gascuel. Maastricht, Pays-Bas: Blackwell Publishers, 1995, pp. 457–468 (cit. on pp. 70, 72).

[Bli82] James F. Blinn. "A Generalization of Algebraic Surface Drawing". In: *ACM Trans. Graph*. 1.3 (July 1982), pp. 235–256 (cit. on p. 70).

[Blo95] Jules Bloomenthal. "Skeletal design of natural forms". PhD thesis. Calgary, Alta., Canada, Canada, 1995 (cit. on p. 58).

[BA11] Craig A Blum and David B Adams. "Who did the first laparoscopic cholecystectomy?" In: *Journal of minimal access surgery* 7.3 (2011), p. 165 (cit. on p. 23).

[CG01] T. Chaperon and F. Goulette. "Extracting cylinders in full 3D data using a random sampling method and the Gaussian image". In: *Vision Modeling and Visualization Conference (VMW'01)*. 2001, pp. 35–42 (cit. on pp. 67, 68).

[CW06] Christian Beder and Wolfgang Förstner. "Direct solutions for computing cylinders from minimal sets of points". In: *9th European Conference on Computer Vision (ECCV'06)*. Vol. 3954. LNCS. May 2006, pp. 135–146 (cit. on p. 67).

[Cle+08] Logan W Clements, William C Chapman, Benoit M Dawant, Robert L Galloway Jr, and Michael I Miga. "Robust surface registration using salient anatomical features for image-guided liver surgery: algorithm and validation". In: *Medical physics* 35.6Part1 (2008), pp. 2528–2540.

[Clo+00] Harry J Cloft, Gregory J Joseph, Frank C Tong, Jonas H Goldstein, and Jacques E Dion. "Use of three-dimensional Guglielmi detachable coils in the treatment of wide-necked cerebral aneurysms". In: *American Journal of Neuroradiology* 21.7 (2000), pp. 1312–1314 (cit. on p. 50).

[CCD10] Olivier Comas, Stéphane Cotin, and Christian Duriez. "A Shell Model for Real-Time Simulation of Intra-ocular Implant Deployment". In: *International Symposium on Computational Models for Biomedical Simulation*. Vol. 5958. Lecture Notes in Computer Science. Springer, 2010, pp. 160–170 (cit. on p. 43).

[Cou+10] Hadrien Courtecuisse, Hoeryong Jung, Jérémie Allard, Christian Duriez, Doo Yong Lee, and Stéphane Cotin. "GPU-based Real-Time Soft Tissue Deformation with Cutting and Haptic Feedback". In: *Progress in Biophysics and Molecular Biology* 103.2-3 (Dec. 2010). Special Issue on Soft Tissue Modelling, pp. 159–168 (cit. on p. 43).
[CHH04] William R Crum, Thomas Hartkens, and DLG Hill. "Non-rigid image registration: theory and practice". In: *The British journal of radiology* 77.suppl_2 (2004), S140-S153 (cit. on p. 125).

[DBB08] Benoit Dagon, Charles Baur, and Vincent Bettschart. "A framework for intra-operative update of 3D deformable models in liver surgery". In: *2008 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society*. IEEE. 2008, pp. 3235-3238 (cit. on p. 118).

[De +11] Hans De Visser, Marcus O Watson, Olivier Salvado, and Joshua D. Passenger. "Progress in virtual reality simulators for surgical training and certification". In: *Medical Journal of Australia* 194.4 (2011), S38 (cit. on p. 2).

[DOG04] Daniel Dekkers, Kees van Overveld, and Rob Golsteijn. "Combining CSG modeling with soft blending using Lipschitz-based implicit surfaces". English. In: *The Visual Computer* 20.6 (2004), pp. 380-391 (cit. on p. 70).

[DNV06] Piet Dirix, Sandra Nuyts, and Walter Van den Bogaert. "Radiation-induced xerostomia in patients with head and neck cancer: a literature review". In: *Cancer: Interdisciplinary International Journal of the American Cancer Society* 107.11 (2006), pp. 2525-2534 (cit. on p. 125).

[Dog13] Issam Doghri. *Mechanics of deformable solids: linear, nonlinear, analytical and computational aspects*. Springer Science & Business Media, 2013 (cit. on p. 26).

[Dou+06] Craig C. Douglas, Jonathan D. Beezley, Janice Coen, Deng Li, Wei Li, Alan K. Mandel, Jan Mandel, Guan Qin, and Anthony Vodacek. "Demonstrating the Validity of a Wildfire DDDAS". In: *Computational Science - ICCS 2006: 6th International Conference, Reading, UK, May 28-31, 2006. Proceedings, Part III*. Springer Berlin Heidelberg, 2006, pp. 522-529.

[Dur+06] C. Duriez, S. Cotin, J. Lenoir, and P. F. Neumann. "New Approaches to Catheter Navigation for Interventional Radiology Simulation". In: *Computer Aided Surgery* 11 (2006), pp. 300-308 (cit. on pp. 30, 31).

[DAK04] Christian Duriez, Claude Andriot, and Abderrahmane Kheddar. "Signorini's contact model for deformable objects in haptic simulations". In: *2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (IEEE Cat. No. 04CH37566)*. Vol. 4. IEEE. 2004, pp. 3232-3237 (cit. on p. 132).

[EP11] Haytham Elhawary and Aleksandra Popovic. "Robust feature tracking on the beating heart for a robotic-guided endoscope". In: *The international journal of medical robotics and computer assisted surgery* 7.4 (2011), pp. 459-468 (cit. on p. 119).

[Fau+12] François Faure, Christian Duriez, Hervé Delingette, Jérémie Allard, Benjamin Gilles, Stéphanie Marchesseau, Hugo Talbot, Hadrien Courtecuisse, Guillaume Bousquet, Igor Peterlik, et al. "Sofa: A multi-model framework for interactive physical simulation". In: *Soft tissue biomechanical modeling for computer assisted surgery*. Springer, 2012, pp. 283-321 (cit. on p. 117).
[Fig+10] Michael Figl, Daniel Rueckert, David Hawkes, Roberto Casula, Mingxing Hu, Ose Pedro, Dong Ping Zhang, Graeme Penney, Fernando Bello, and Philip Edwards. "Image guidance for robotic minimally invasive coronary artery bypass". In: *Computerized Medical Imaging and Graphics* 34.1 (2010), pp. 61–68 (cit. on p. 114).

[Fis+02] Bruce Fischl, David H Salat, Evelina Busa, Marilyn Albert, Megan Dieterich, Christian Haselgrove, Andre Van Der Kouwe, Ron Killiany, David Kennedy, Shuna Klaveness, et al. "Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain". In: *Neuron* 33.3 (2002), pp. 341–355 (cit. on p. 54).

[FU06] J. Florens and D. Urma. "Dynamical Issues at the Low Level of Human/Virtual Object Interaction". In: *2006 14th Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems*. Mar. 2006, pp. 315–320 (cit. on p. 26).

[Fri+10] O. Friman, M. Hindennach, C. Kühnel, and H.-O. Peitgen. "Multiple hypothesis template tracking of small 3D vessel structures". In: *MedIA 14.2 (2010)*, pp. 160–171 (cit. on pp. 67, 77).

[Gol05] Ron Goldman. "Curvature formulas for implicit curves and surfaces". In: *Comput. Aided Geom. Des.* 22.7 (Oct. 2005), pp. 632–658 (cit. on p. 72).

[HL00] E Mark Haacke and Zhi-Pei Liang. "Challenges of imaging structure and function with MRI". In: *IEEE Engineering in Medicine and Biology Magazine* 19.5 (2000), pp. 55–62 (cit. on p. 54).

[HK76] J.C. Halpin and J.L. Kardos. "The halpin-tsai equations: a review". In: *Polymer Engineering Science* 16 (1976), pp. 344–52 (cit. on p. 33).

[Hao15] Nazim Haouchine. "Image-guided simulation for augmented reality during hepatic surgery". PhD thesis. Université des Sciences et Technologie de Lille - Lille I, 2015 (cit. on pp. 5, 123).

[HK98] SA Hojjatoleslami and Josef Kittler. "Region growing: a new approach". In: *IEEE Transactions on Image processing* 7.7 (1998), pp. 1079–1084 (cit. on p. 54).

[Hum03] Jay D Humphrey. "Review Paper: Continuum biomechanics of soft biological tissues". In: *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*. Vol. 459. The Royal Society. 2003, pp. 3–46 (cit. on p. 2).

[Iwa+11] Naoya Iwamoto, Ryusuke Sagawa, Shoji Kunitomo, and Shigeo Morishima. "Estimating fluid simulation parameters from videos". In: *ACM SIGGRAPH 2011 Posters*. 2011, p. 3.

[Jel+02] James S Jelinek, Mark D Murphey, James A Welker, Robert M Henshaw, Mark J Kransdorf, Barry M Shmookler, and Martin M Malawer. "Diagnosis of Primary Bone Tumors with Image-guided Percutaneous Biopsy: Experience with 110 Tumors 1". In: *Radiology* 223.3 (2002), pp. 731–737 (cit. on p. 1).
[Jia17] Zhifan Jiang. "Évaluation des mobilités et modélisation géométrique du système pelvien féminin par analyse d'images médicales". PhD thesis. Université des Sciences et Technologie de Lille, 2017 (cit. on pp. 5, 88).

[KWT88] Michael Kass, Andrew Witkin, and Demetri Terzopoulos. "Snakes: Active contour models". In: *International journal of computer vision* 1.4 (1988), pp. 321–331 (cit. on p. 54).

[Ker06] Amy E Kerdok. "Characterizing the nonlinear mechanical response of liver to surgical manipulation". In: (2006) (cit. on p. 115).

[Kim+12] Jae-Hak Kim, Adrien Bartoli, Toby Collins, and Richard Hartley. "Tracking by detection for interactive image augmentation in laparoscopy". In: *International Workshop on Biomedical Image Registration*. Springer. 2012, pp. 246–255 (cit. on p. 114).

[KP93] S Kumar and A Petho. "An algorithm for the numerical inversion of a tridiagonal matrix". In: *Communications in Numerical Methods in Engineering* 9.4 (1993), pp. 353–359 (cit. on p. 33).

[Lec+15] Pauline Lecomte-Grosbras, J-F Witz, Mathias Brieu, N Faye, Michel Cosson, and Chrystele Rubod. "Quantification of pelvic mobility on dynamic magnetic resonance images: using mechanical insight to help diagnose pelvic pathologies". In: *Strain* 51.4 (2015), pp. 301–310 (cit. on p. 108).

[Len+05] Julien Lenoir, Stéphane Cotin, Christian Duriez, and Paul Neumann. "Interactive physically-based simulation of Catheter and Guidewire". In: *Second Workshop in Virtual Reality Interactions and Physical Simulations (VRIPHYS)*. Pisa - Italy, Nov. 2005 (cit. on p. 33).

[LK+81] Bruce D Lucas, Takeo Kanade, et al. "An iterative image registration technique with an application to stereo vision". In: (1981) (cit. on pp. 101, 102, 119).

[Mar+10] Stéphanie Marchesseau, Tobias Heimann, Simon Chatelin, Rémy Willinger, and Hervé Delingette. "Fast porous visco-hyperelastic soft tissue model for surgery simulation: application to liver surgery". In: *Progress in biophysics and molecular biology* 103.2-3 (2010), pp. 185–196 (cit. on p. 116).

[Mar+97] J. A. Martin, G. Regehr, R. Reznick, H. MacRae, J. Murnaghan, Hutchison C., and Brown M. "Objective Structured Assessment of Cataract Surgical Skill". In: *The British journal of surgery* 84 (2 1997), pp. 273–278 (cit. on p. 46).

[Mor+09] Francesc Moreno-Noguer, Mathieu Salzmann, Vincent Lepetit, and Pascal Fua. "Capturing 3D stretchable surfaces from single images in closed form". In: *2009 IEEE Conference on Computer Vision and Pattern Recognition*. IEEE. 2009, pp. 1842–1849 (cit. on p. 119).

[Mul+10] Patrick Mullen, Fernando De Goes, Mathieu Desbrun, David Cohen-Steiner, and Pierre Alliez. "Signing the unsigned: Robust surface reconstruction from raw pointsets". In: *Computer Graphics Forum*. Vol. 29. 5. Wiley Online Library. 2010, pp. 1733–1741 (cit. on p. 63).
[Mur91] S. Muraki. "Volumetric shape description of range data using Blobby Model". In: *SIGGRAPH Comput. Graph.* 25 (4 1991), pp. 227-235 (cit. on pp. 70, 72, 73).

[Nea+06] Andrew Nealen, Matthias Müller, Richard Keiser, Eddy Boxerman, and Mark Carlson. "Physically based deformable models in computer graphics". In: *Computer graphics forum*. Vol. 25. 4. Wiley Online Library. 2006, pp. 809-836 (cit. on p. 26).

[Nel+12] Benjamin E Nelms, Greg Robinson, Jay Markham, Kyle Velasco, Steve Boyd, Sharath Narayan, James Wheeler, and Mark L Sobczak. "Variation in external beam treatment plan quality: an inter-institutional study of planners and planning systems". In: *Practical radiation oncology* 2.4 (2012), pp. 296-305 (cit. on p. 125).

[NPF05] Matthieu Nesme, Yohan Payan, and François Faure. "Efficient, physically plausible finite elements". In: 2005 (cit. on p. 116).

[Nic+11] Stéphane Nicolau, Luc Soler, Didier Mutter, and Jacques Marescaux. "Augmented reality in laparoscopic surgical oncology". In: *Surgical Oncology* 20.3 (2011), pp. 189-201 (cit. on p. 113).

[Oht+05] Yutaka Ohtake, Alexander Belyaev, Marc Alexa, Greg Turk, and Hans-Peter Seidel. "Multi-level partition of unity implicits". In: *Acm Siggraph 2005 Courses*. 2005, 173-es (cit. on p. 63).

[Okt+13] Ozan Oktay, Li Zhang, Tommaso Mansi, Peter Mountney, Philip Mewes, Stéphane Nicolau, Luc Soler, and Christophe Chefd'hotel. "Biomechanically driven registration of pre-to intra-operative 3d images for laparoscopic surgery". In: *International Conference on Medical Image Computing and Computer-Assisted Intervention*. Springer. 2013, pp. 1-9 (cit. on p. 118).

[PDC11] Igor Peterlik, Christian Duriez, and Stéphane Cotin. "Asynchronous haptic simulation of contacting deformable objects with variable stiffness". In: *2011 IEEE/RSJ International Conference on Intelligent Robots and Systems*. IEEE. 2011, pp. 2608-2613 (cit. on p. 37).

[PDC12] Igor Peterlik, Christian Duriez, and Stéphane Cotin. "Modeling and real-time simulation of a vascularized liver tissue". In: *International Conference on Medical Image Computing and Computer-Assisted Intervention*. Springer. 2012, pp. 50-57 (cit. on pp. 116, 117).

[Pha01] Dzung L Pham. "Robust fuzzy segmentation of magnetic resonance images". In: *Computer-Based Medical Systems, 2001. CBMS 2001. Proceedings*. 14th IEEE Symposium on. IEEE. 2001, pp. 127-131 (cit. on p. 54).

[Pra+10] Philip Pratt, Danail Stoyanov, Marco Visentini-Scarzanella, and Guang-Zhong Yang. "Dynamic guidance for robotic surgery using image-constrained biomechanical models". In: *International Conference on Medical Image Computing and Computer-Assisted Intervention*. Springer. 2010, pp. 77-85 (cit. on p. 114).

[Prz68] J.S. Przemieniecki. *Theory of Matrix Structural Analysis*. McGraw-Hill, 1968 (cit. on pp. 30, 116).
[PM13] GA Puerto-Souza and GL Mariottini. "An augmented-reality system for laparoscopic surgery robust to complete occlusions and fast camera motions". In: *International Conference on Robotics and Automation*. 2013 (cit. on p. 114).

[RF14] Julien Réthoré and Marc François. "Curve and boundaries measurement using B-splines and virtual images". In: *Optics and Lasers in Engineering* 52 (2014), pp. 145-155 (cit. on p. 91).

[RKK03] Josef Rösch, Frederick S Keller, and John A Kaufman. "The birth, early years, and future of interventional radiology". In: *Journal of Vascular and Interventional Radiology* 14.7 (2003), pp. 841-853 (cit. on p. 22).

[Rot82] Scott D Roth. "Ray casting for modeling solids". In: *Computer graphics and image processing* 18.2 (1982), pp. 109-144 (cit. on p. 120).

[Sal+07] Georges M. Saleh, V. Gauba, Arijit Mitra, Andre S. Litwin, Andrew K. K. Chung, and Larry Benjamin. "Objective Structured Assessment of Cataract Surgical Skill". In: *Archives of Ophthalmology* 125 (3 2007), pp. 363-366 (cit. on p. 46).

[SF10] Mathieu Salzmann and Pascal Fua. "Linear local models for monocular reconstruction of deformable surfaces". In: *IEEE Transactions on Pattern Analysis and Machine Intelligence* 33.5 (2010), pp. 931-944 (cit. on p. 119).

[Sat+98] K Satoh, J Satomi, S Matsubara, and S Nagahiro. "Measurement of volume ratio to predict coil compaction, on aneurysmal embolization". In: *Interventional Neuroradiology* 4.1_suppl (1998), pp. 179-182 (cit. on p. 50).

[Sch22] Pierre Schegg. "Autonomous Catheter and Guidewire Navigation for Robotic Percutaneous Coronary Interventions". PhD thesis. Université des Sciences et Technologie de Lille, 2022 (cit. on p. 5).

[SFA10] Benoit Semin, Marc Louis Maurice François, and Harold Auradou. "Analytical shape determination of fiber-like objects with Virtual Image Correlation". In: *arXiv preprint arXiv:1001.2636* (2010) (cit. on p. 91).

[SV16] N Senthilkumaran and S Vaithegi. "Image segmentation by using thresholding techniques for medical images". In: *Computer Science & Engineering: An International Journal* 6.1 (2016), pp. 1-13 (cit. on p. 54).

[Ser+03] Maxime Sermesant, Clément Forest, Xavier Pennec, Hervé Delingette, and Nicholas Ayache. "Deformable biomechanical models: Application to 4D cardiac image analysis". In: *Medical image analysis* 7.4 (2003), pp. 475-488 (cit. on p. 114).

[SLH10] Tian Shen, Hongsheng Li, and Xiaolei Huang. "Active volume models for medical image segmentation". In: *IEEE transactions on medical imaging* 30.3 (2010), pp. 774-791 (cit. on p. 114).

[She99] A. Sherstyuk. "Kernel Functions in Convolution Surfaces: A Comparative Analysis". In: *The Visual Computer* 15.4 (1999), pp. 171-182 (cit. on p. 70).

[Sig33] Antonio Signorini. "Sopra alcune questioni di statica dei sistemi continui". In: *Annali della Scuola Normale Superiore di Pisa-Classe di Scienze* 2.2 (1933), pp. 231-251 (cit. on p. 132).
[Sil+09] M Silvennoinen, J-P Mecklin, P Saariluoma, and T Antikainen. "Expertise and skill in minimally invasive surgery". In: *Scandinavian Journal of Surgery* 98.4 (2009), pp. 209-213 (cit. on p. 1).

[Sin+02] V. Singha, D. Gressa, R. Higashidab, C. Dowdb, V. Halbachb, and S. Johnston. "The Learning Curve for Coil Embolization of Unruptured Intracranial Aneurysms". In: *American Journal of Neuroradiology* 23 (2002), pp. 768-771 (cit. on p. 47).

[Spe+11] S Speidel, S Roehl, S Suwelack, R Dillmann, H Kenngott, and B Mueller-Stich. "Intraoperative surface reconstruction and biomechanical modeling for soft tissue registration". In: *Proc. Joint Workshop on New Technologies for Computer/Robot Assisted Surgery*. 2011, pp. 0-0 (cit. on p. 114).

[Ste13] Martin Oliver Steinhauser. *Computer simulation in physics and engineering*. Walter de Gruyter, 2013 (cit. on p. 21).

[Su+09] Li-Ming Su, Balazs P Vagvolgyi, Rahul Agarwal, Carol E Reiley, Russell H Taylor, and Gregory D Hager. "Augmented reality during robot-assisted laparoscopic partial nephrectomy: toward real-time 3D-CT to stereoscopic video registration". In: *Urology* 73.4 (2009), pp. 896-900 (cit. on p. 114).

[Tal14] Hugo Talbot. "Interactive Patient-Specific Simulation of Cardiac Electrophysiology". PhD thesis. Université des Sciences et Technologies de Lille, 2014.

[Tam+02] Shinichi Tamatani, Yasushi Ito, Hiroshi Abe, Tetsuo Koike, Shigekazu Takeuchi, and Ryuichi Tanaka. "Evaluation of the stability of aneurysms after embolization using detachable coils: correlation between stability of aneurysms and embolized volume of aneurysms". In: *American Journal of Neuroradiology* 23.5 (2002), pp. 762-767 (cit. on p. 50).

[Tau91] G. Taubin. "Estimation of Planar Curves, Surfaces, and Nonplanar Space Curves Defined by Implicit Equations with Applications to Edge and Range Image Segmentation". In: *PAMI* 13 (11 1991), pp. 1115-1138 (cit. on pp. 73, 81, 82).

[Teb+09] Dogu Teber, Selcuk Guven, Tobias Simpfendörfer, Mathias Baumhauer, Esref Oguz Güven, Faruk Yencilek, Ali Serdar Gözen, and Jens Rassweiler. "Augmented reality: a new tool to improve surgical accuracy during laparoscopic partial nephrectomy? Preliminary in vitro and in vivo results". In: *European urology* 56.2 (2009), pp. 332-338 (cit. on p. 114).

[Fel00] C. A. Felippa. *A systematic approach to the element independent corotational dynamics of finite elements*. Tech. rep. CU-CAS-00-03. Center for Aerospace Structures, 2000 (cit. on pp. 31, 116).

[TBC95] N. Tsingos, E. Bittar, and M.-P. Cani. "Implicit Surfaces for Semi-Automatic Medical Organ Reconstruction". In: *Computer Graphics Internat. (CGI'95)*. 1995, pp. 3-15 (cit. on pp. 72, 73).
[Vei+14] Catarina Veiga, Jamie McClelland, Syed Moinuddin, Ana Lourenço, Kate Rick- etts, James Annkah, Marc Modat, Sébastien Ourselin, Derek D'Souza, and Gary Royle. "Toward adaptive radiotherapy for head and neck patients: feasibility study on using CT-to-CBCT deformable registration for "dose of the day" calculations". In: *Medical physics* 41.3 (2014) (cit. on p. 125).

[WMK06] Ruikang K Wang, Zhenhe Ma, and Sean J Kirkpatrick. "Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue". In: *Applied Physics Letters* 89.14 (2006), p. 144103 (cit. on p. 21).

[Web03] Steve Webb. "The physical basis of IMRT and inverse planning". In: *The British journal of radiology* 76.910 (2003), pp. 678-689 (cit. on p. 125).

[WNV00] Onno Wink, Wiro J. Niessen, and Max A. Viergever. "Fast Delineation and Visualization of Vessels in 3D Angiographic Images". In: *IEEE Trans. Med. Imaging* 19.4 (2000), pp. 337-346 (cit. on p. 69).

[WH94] Andrew P Witkin and Paul S Heckbert. "Using particles to sample and control implicit surfaces". In: *Proceedings of the 21st annual conference on Computer graphics and interactive techniques*. 1994, pp. 269-277.

[WC07a] W.C.K Wong and A.C.S Chung. "Probabilistic vessel axis tracing and its application to vessel segmentation with stream surfaces and minimum cost paths". In: *MedIA* 11 (2007), pp. 567-587 (cit. on pp. 67, 68).

[WC07b] Wilbur CK Wong and Albert CS Chung. "Probabilistic vessel axis tracing and its application to vessel segmentation with stream surfaces and minimum cost paths". In: *Medical Image Analysis* 11.6 (2007), pp. 567-587.

[WLS12] Stefanie Wuhrer, Jochen Lang, and Chang Shu. "Tracking complete deformable objects with finite elements". In: *2012 Second International Conference on 3D Imaging, Modeling, Processing, Visualization & Transmission*. IEEE. 2012, pp. 1-8 (cit. on p. 119).

[YAG+05] Kenji YAGI, Koichi SATOH, Junichiro SATOMI, Shunji MATSUBARA, and Shinji NAGAHIRO. "Evaluation of Aneurysm Stability After Endovascular Embolization With Guglielmi Detachable Coils: Correlation Between Long-Term Stability and Volume Embolization Ratio". In: *Neurologia medico-chirurgica* 45.11 (2005), pp. 561-566 (cit. on p. 50).

[YC14] Michael C Yip and David B Camarillo. "Model-less feedback control of continuum manipulators in constrained environments". In: *IEEE Transactions on Robotics* 30.4 (2014), pp. 880-889 (cit. on p. 133).

[Yur14] Ahmed Yureidini. "Robust blood vessel surface reconstruction for interactive simulations from patient data". PhD thesis. Université des Sciences et Technologie de Lille - Lille I, 2014 (cit. on p. 5).

[Zac01] Gabriel Zachmann. "Optimizing the collision detection pipeline". In: *Proc. of the First International Game Technology Conference (GTEC)*. Vol. 2. 2001 (cit. on p. 56).
[Zha19] Zhongkai Zhang. "vision-based calibration, position control and force sensing for soft robots". PhD thesis. université de lille, Jan. 2019 (cit. on p. 5).

[ZZ87] Olgierd C Zienkiewicz and Jian Z Zhu. "A simple error estimator and adaptive procedure for practical engineering analysis". In: *International journal for numerical methods in engineering* 24.2 (1987), pp. 337–357 (cit. on p. 34).
List of Figures

<table><tr><td>1.1</td><td>Photography of a radiology intervention: thin flexible devices are inserted in the femoral artery and are navigated inside the vascular network until reaching the diseased area. Then different treatments can be performed. Monitoring the motion of the device in the patient is performed via X-ray imaging. Photography appears courtesy of www.dicardiology.com</td><td>23</td></tr><tr><td>1.2</td><td>Illustrations of laparoscopic surgery procedure: 4 small incisions are required to insert surgical instruments and camera (left). The surgeons then monitor the procedure on a video screen while manipulating the long thin surgical instruments (right). Left image appears courtesy of Sunshine Coast Advanced Laparoscopic Surgery, right image appears courtesy of www.chronogram.com.</td><td>24</td></tr><tr><td>1.3</td><td>Example of a neuro-surgical training system by Neurotouch©. The system includes robotic arm equipped with stylus (Geomagic© Touch) to acquire user gesture and to render forces in case of contact and a stereo-vision display to mimic the use of binoculars.</td><td>25</td></tr><tr><td>1.4</td><td>A simplified overview of the simulation loop for bio-mechanical applications. These 4 steps are repeated to produce feedbacks (visual, haptic, audio...) to the user.</td><td>26</td></tr><tr><td>1.5</td><td>Examples of manufactured surgical devices than can be modeled as one-dimensional solids: Surgical Thread (®AD Surgical), Endoscope (®International Light Technologies), Catheter (®Paralogic), Detachable coil (®MedTronic)</td><td>28</td></tr><tr><td>1.6</td><td>Example of a coiling procedure in an aneurysm (blood vessel bulge): a thin wire (coil) is inserted in the aneurysm reducing the blood flow thus the blood pressure and allowing the blood to clot.</td><td>29</td></tr><tr><td>1.7</td><td>Examples of the diversity of manufactured detachable coils used in interventional radiology. Length, diameter and shapes (helical, 3D, 360 degrees) may differ and are used given the pathology encountered.</td><td>29</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>1.8</td>
      <td>Aneurysm coiling requires several coils to fill the aneurysm in order to reduce the blood flow and thus prevents the aneurysm from rupturing. Image courtesy of Stryker Neurovascular.</td>
      <td>30</td>
    </tr>
    <tr>
      <td>1.9</td>
      <td>Modeling of a wire using linked beam elements. (a) Given an initial rest shape and forces applied on DOFs, internal forces are computed using a local reference frame to handle large deformations. Starting from the first beam, reference frames are updated to suppress rigid transformation from the expression of the beam deformation. (b) Modification of the initial shape by moving frames or positions allows the generation of various deformed shapes. (c) A composite model (co-axial combination of two wires) is obtained by geometrically mixing different rest shapes on one geometric model.</td>
      <td>32</td>
    </tr>
    <tr>
      <td>1.10</td>
      <td>Experimental setup: (a) the coil is deployed in a contact-free environment. (b) Volumetric data is obtained by 3DXA and a marching cube is performed to get a mesh of the real coil. (c) Central line of the coil is extracted and a continuous formulation is built using splines.</td>
      <td>34</td>
    </tr>
    <tr>
      <td>1.11</td>
      <td>Visual comparison of our coil simulation (in red) with the reconstructed coil model from our experiments (in yellow) at different stages of the deployment. Quantitative errors are available in table 1.1.</td>
      <td>35</td>
    </tr>
    <tr>
      <td>1.12</td>
      <td>Asynchronous simulation: details on the steps performed by the two loops that run at different frequencies. The central blocks represent shared buffers that are needed to compute mechanical interactions.</td>
      <td>38</td>
    </tr>
    <tr>
      <td>1.13</td>
      <td>Simulation of a deformable cylinder (simulated at 50 Hz) and a thread (simulated at 1 kHz) interacting with a solid plane.</td>
      <td>39</td>
    </tr>
    <tr>
      <td>1.14</td>
      <td>Functional view and picture of the training system. The system relies on devices used in real procedures (microscope, pedals...), replicas of surgical tools that are tracked in real-time and on a simulation engine that renders a realistic and interactive simulation of cataract surgery.</td>
      <td>40</td>
    </tr>
    <tr>
      <td>1.15</td>
      <td>Tools replicas</td>
      <td>41</td>
    </tr>
    <tr>
      <td>1.16</td>
      <td>Overview of our setup: mannequin and microscope where 6 infrared cameras are attached. These cameras allow to track in real-time the motion of the surgical tools.</td>
      <td>42</td>
    </tr>
    <tr>
      <td>1.17</td>
      <td>Screenshots of the tracking of surgical tools. With the reflective markers, the tracking system is able to compute the 3D position of the tool in real-time and using different marker signatures, we can manage to detect several tools simultaneously without ambiguity.</td>
      <td>42</td>
    </tr>
    <tr>
      <td>1.18</td>
      <td>Tessellation algorithm that can dynamically refine tetrahedra (two images on the left). Two-levels tetrahedral topology for lens phacoemulsification (two pictures on the right).</td>
      <td>43</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>1.19</td>
      <td>Lens implant and its mesh</td>
      <td>44</td>
    </tr>
    <tr>
      <td>1.20</td>
      <td>Three screenshots of the phacoemulsification simulation. Local fragmentation of the lens create a cross-like pattern on the lens.</td>
      <td>45</td>
    </tr>
    <tr>
      <td>1.21</td>
      <td>From left to right: tearing propagation throughout the tissue. Fibers (not visible on the meshes since they are defined per triangles) are oriented with a 45 degree angle, and forces are applied outward, on the left and right sides of the square mesh. Our fracture criteria is represented as the shaded area. One can see on the left image how the anisotropy affects the deformation of the tissue, and on the rightmost images how the fracture direction depends on the previous fracture direction, as well as fiber direction. The square mesh is composed of 1,500 triangles. Anisotropic deformation, fracture and remeshing are performed in real-time.</td>
      <td>46</td>
    </tr>
    <tr>
      <td>1.22</td>
      <td>Simulation of capsulorhexis, a technique used to create a circular opening in the lens capsule. This technique relies essentially on the application of shear and stretch forces to propagate a fracture throughout the membrane. Top (a): the mesh used to support the computation of the deformation and topological changes. Bottom (b): final result as seen in the training system.</td>
      <td>47</td>
    </tr>
    <tr>
      <td>1.23</td>
      <td>Reconstructed vasculature using 3DXA imaging modality. The 3D volume was segmented and reconstructed using triangles.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>1.24</td>
      <td>Examples of our simulation results: (two images on the left) real coil embolization; (two images on the right) our simulated coil embolization with 3D coils.</td>
      <td>49</td>
    </tr>
    <tr>
      <td>1.25</td>
      <td>Examples of our simulation results: close-up on the coil deployment and wire-frame display to illustrate the deformation of the coil.</td>
      <td>49</td>
    </tr>
    <tr>
      <td>1.26</td>
      <td>Examples of our simulation results: pressure of the coil on the aneurysm surface can be displayed using a thermal scale Using a larger coil (right) induces more important pressure on the aneurysm wall which could lead to rupture.</td>
      <td>50</td>
    </tr>
    <tr>
      <td>1.27</td>
      <td>Examples of our simulation results: using helical coils with different diameters: (left) 4mm helical coil which is too small for the aneurysm and the coil is floating inside the aneurysm and this can be dangerous for the patient, (middle) 7mm helical coil which is the optimal size for this aneurysm, (right) a 10mm helical coil which is too big and leads to extra pressure on the aneurysm wall.</td>
      <td>51</td>
    </tr>
    <tr>
      <td>1.28</td>
      <td>Computation of the Coil Embolization Surface Ratio metric: it divides the surface of the coil by the surface of the aneurysm.</td>
      <td>52</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>2.1</td>
      <td>200msec recorded sequences of catheter deployments in a small cylindrical vessels. Catheters share the same mechanical parameters and only the friction coefficient between the vessel wall and the tool is modified. This illustrates the importance of adequate friction estimation and smooth geometry.</td>
      <td>55</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Recorded sequence of a coil deployed at low and constant speed where image acquisition is performed at 15 fps through 3D Rotational Angiography. The sequence shows the effect of friction and then sliding as the coil completely changes its configuration between frames 3 and 4 (67 msecs between each frame).</td>
      <td>56</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Reconstructed liver and its vasculature: segmentation and reconstruction have been performed without smoothing. The limited resolution of the imaging device leads to virtual 3D models that exhibit jagged surfaces which are not very organic.</td>
      <td>57</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Skeleton-based convolution surface: using a skeleton defined by cylindrical-shaped primitives (top, for clarity, only the center line of the cylinders are shown) and a Gaussian convolution filter. The potential field is displayed with a heat-map (left-bottom).</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Example of catheter and guide-wire navigation inside a vascular phantom using convolution surfaces as the underlying geometric model.</td>
      <td>60</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Example of catheter and guide-wire navigation inside a vascular network of a specific patient. The vascular model consists in more than 4000 blood vessels and the simulation can be run at interactive rates (100 Hz). Complex configurations (tortuous and small vessels) are handled correctly.</td>
      <td>61</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Catheter deployment in a closed capsule. (Left-column) the capsule is modeled as a polygonal mesh, (Right-column) the capsule is modeled as a convolution surface. The snapshots taken at the same simulation time-steps illustrate that artificial friction can be induced with polygonal meshes because the tip sticks into the perfectly smooth vessel instead of sliding along it (which is the case of the implicit surface).</td>
      <td>62</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Examples of reconstruction issues from segmented medical images. Since the contouring of the vascular structure do not encode any topological information many artifacts may occur such as blending between aneurysms and vessels (top) or holes and discontinuing vessels (bottom).</td>
      <td>64</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>2.9</td>
      <td>Outline of RBT vessel tracking algorithm. (top left) Cut plane through the original 3DRA data; (top right) Prediction: Starting from C<sub>t-1</sub> (green dot) and d&#773;<sub>t-1</sub> (white arrow), the new center C<sub>t</sub><sup>*</sup> (red dot) is found. Candidate points (white dots) are extracted at the vessel surface by casting rays in the volume, from C<sub>t</sub><sup>*</sup>. Only the points around the displayed cut plane are shown. (bottom left) N<sub>d</sub> directions are tested for the cylinder axes (yellow arrows). One fitting cylinder is found per axis, using RANSAC. (bottom right) Best cylinder found (C<sub>t</sub>, r<sub>t</sub>, d&#773;<sub>t</sub>), together with its consensus set (green dots). The red dots were used to compute the cylinder height h<sub>t</sub>: they are the set of points P<sub>t</sub>, attached to the cylinder.</td>
      <td>66</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>Adding a new branch: (left) Initialize with a user-input point (green point), project it on the parent centerline (violet point). Their connection initializes the vessel direction (arrow). Tune the initial (green) point, by moving it along the vessel direction, so that it is out of the parent vessel, but still close to the proximal end of the child vessel (blue point). (right) Tracks the child branch from this point and connects it to the parent branch (stale blue polyline).</td>
      <td>69</td>
    </tr>
    <tr>
      <td>2.11</td>
      <td>Outline of the LIM algorithm. (top left) Input data: RBT for two branching arteries on 3DRA data results in points on the centerlines (blue and black points) equipped with points locally extracted at the vessel surface (orange dots); (top right) The red node is chosen. (middle left) Closest neighbor nodes are selected w.r.t a given topological distance (of one, here). (middle right) Local point sets are concatenated. (bottom left) The local implicit surface (blobby model) is initialized with two close spheres. (bottom right) LIM successively fits and splits blobs to fit the input points. 30 subdivisions were necessary to produce the resulting implicit surface.</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.12</td>
      <td>The elementary reconstruction step in 2D. (a) Points (black) are to be reconstructed with two blobs (red contour). (b) First, all blobs are fine tuned over the centers and then, over the widths. (c) Second, find the point farthest to the contour (cyan) and the blob closest to it (purple). (d) Third, replace this blob with two blobs. The first (black contour) is placed at the same position (black) and the second (cyan contour) is translated in the direction of the farthest point (cyan). Widths are calculated so that both blobs centered at the black point generate the same isosurface as the fissioned blob (purple). (e) Fourth, the translated blob (cyan) is tuned over its center then width. (f) Finally, the BM is updated by integrating the recently tuned blob.</td>
      <td>74</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>2.13</td>
      <td>Implicit modeling of an aneurysm. The points {P<sub>i</sub>} are in red. (From left to right) Initialization with a single blob; after the first minimization; after 25 subdivisions; final result (100 subdivisions)</td>
      <td>76</td>
    </tr>
    <tr>
      <td>2.14</td>
      <td>Top: (left) Maximum Intensity Projection of one patient data. Center-line delineation produced by MHT (middle) and RBT (right). Bottom: (left) When the MHT centerline (red) does not remain inside the vessel wall, it is cut – for TL computation – before the problem arises (green). (middle) MHT centerline (red) and the points (cyan) selected for MASD computation (w = r̄). (right) The resulting points for MHT (red) and RBT (blue, all the original points) used for MASD computation.</td>
      <td>79</td>
    </tr>
    <tr>
      <td>2.15</td>
      <td>(left) A tortuous branch of the posterior choroidal artery (0.6 – 0.7 mm of diameter) in Patient 3. (middle) MHT is misled by neighboring structures in the proximal portion of the artery whereas (right) RBT tracks the tortuous vessel and successfully handles the Kissing Vessel (KV) issue.</td>
      <td>79</td>
    </tr>
    <tr>
      <td>2.16</td>
      <td>Reference surface (blue) and the resulting BM (purple) for maximal values in Table 2.2. The 300 noisy input points are in red.</td>
      <td>81</td>
    </tr>
    <tr>
      <td>2.17</td>
      <td>Visual assessment of LIM: (top-left) isosurface of 3DRA data; (top-right) LIM result, showing much more vessels; (bottom-left) no unwanted blending at kissing vessels (encircled); (bottom-right) difficulty to sometimes model branching of a tiny vessel on a large one (encircled), small discontinuities are visible between adjacent BMs (arrow).</td>
      <td>83</td>
    </tr>
    <tr>
      <td>2.18</td>
      <td>Data used for validation: capsule (left) and spiral (right) surfaces. Both synthetic shapes presented respectively dimensions (width×height×depth) in mm: 11.9 × 2.5 × 2.5 and 8.1 × 6.12 × 1.05. A primitive was placed at each red dot composing the skeleton (white).</td>
      <td>85</td>
    </tr>
    <tr>
      <td>2.19</td>
      <td>The surface of the local BM selected to solve contact at the tool tip point is displayed in wire-frame for 4 simulation steps. The overall vessel surface is shown in transparent red as a mesh.</td>
      <td>87</td>
    </tr>
    <tr>
      <td>2.20</td>
      <td>Sagittal dynamic magnetic resonance (MR) image (spatial resolution = 1.17 mm/pixel, 256*256 pixels. (top): First and last images of a dynamic MRI sequence of a healthy volunteer. (bottom-left): dynamic MRI of a prolapsus where the uterus moved down. (bottom-right): static MRI of endometriosis</td>
      <td>90</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>2.21</td>
      <td>Contour detection with the cost functions. (A) Three small circles initialised by user with 10 control points for each. (B) Positions after an affine transformation. (C) Final contours with optimal number of control points. (D-F) Illustration of the cost function metrics E<sub>image</sub> along the organ contours: [I<sub>r</sub>(u, v) - I<sub>V</sub>(u, v)]<sup>2</sup>, with u evolves in the tangent direction of the contour and v in the direction perpendicular to the contour.</td>
      <td>95</td>
    </tr>
    <tr>
      <td>2.22</td>
      <td>Measurement of the detection error between a correct contour (green) and an arbitrary contour (red): (A) Hausdorff distance = length of s. (B) Average curve distance = average length of s. (C) Dice coefficient = s/(s + t).</td>
      <td>98</td>
    </tr>
    <tr>
      <td>2.23</td>
      <td>Illustration of results: green curve shows the manually segmented organ as the reference, red curve shows the result of semi-automatic algorithm and yellow line shows the Hausdorff distance. (A) Witness patient (1.17 mm/pixel). (B) Patient with endometriosis (0.49 mm/pixel). (C) Patient with prolapse (1.17 mm/pixel).</td>
      <td>98</td>
    </tr>
    <tr>
      <td>2.24</td>
      <td>Different exceptions of detection caused by unwanted effect in images: (A) As rectum may be folded somewhere, wrinkles appear in image as black rings. (B) The injected gel and intrinsic material may produce blur parts in image. (C) Intrinsic material show spots in rectum.</td>
      <td>99</td>
    </tr>
    <tr>
      <td>2.25</td>
      <td>(left-middle) Choice of the features along the B-spline contour. (right) Adjusted model at instant t<sub>0</sub>, The figure is generated by our application.</td>
      <td>102</td>
    </tr>
    <tr>
      <td>2.26</td>
      <td>Tracking of the features using the optical flow method: window grid size = 30 × 30 pixel.</td>
      <td>103</td>
    </tr>
    <tr>
      <td>2.27</td>
      <td>Deformation of the model from instant t to t + 1: firstly, the features are deformed using optical flow tracking method, for example the transformations of point 0 and 1 are computed by the optical flow equation. Then, the same transformations are applied on their neighboring free control points to form the intermediate model (instant t*). Finally, the free control points (for example the ones situated between 0 et 1) are adjusted in the direction perpendicular to the segment using virtual image correlation, and produce the model at instant t + 1.</td>
      <td>104</td>
    </tr>
    <tr>
      <td>2.28</td>
      <td>Computation of the displacements (A) and shears (B), between instant t<sub>0</sub> and t.</td>
      <td>105</td>
    </tr>
    <tr>
      <td>2.29</td>
      <td>Comparison of the angles of rotation computed by two methods: the results given by Bézier method are shown on the top and the ones given by B-spline method are on the bottom. The white contour is the initial one. The imposed rotations are 2.5, 10.0 and 17.5 degrees. The color map illustrate the variation along the contour in two cases.</td>
      <td>107</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>2.30</td>
      <td>Comparison of the angles of rotation computed by two methods: the results given by Bézier method are shown on the top and the ones given by B-spline method are on the bottom. The white contour is the initial one. The imposed rotations are 2.5, 10.0 and 17.5 degrees. The color map illustrate the variation along the contour in two cases.</td>
      <td>108</td>
    </tr>
    <tr>
      <td>2.31</td>
      <td>Illustration of visual comparison between two method. On the top: displacements at the middle of straining computed by (A) Bézier method and (B) B-spline method. The white contours are the initial ones. On the bottom: (C) initial position of segments illustrating the absence of shear, and deviation of the segments at the middle of straining computed by (D) Bézier method and (E) B-spline method.</td>
      <td>109</td>
    </tr>
    <tr>
      <td>2.32</td>
      <td>Quantitative comparison of the displacements computed by two methods: (A,B) Displacement computed by B-spline method at the middle and end of straining, (C,D) Displacement computed by Bézier method at the middle and end of straining.</td>
      <td>110</td>
    </tr>
    <tr>
      <td>2.33</td>
      <td>Quantitative evaluation of the displacement field in the context of hyper-mobility: (A) Displacements at the beginning of straining and (B) at the middle of straining. The white contours are the initial ones.</td>
      <td>111</td>
    </tr>
    <tr>
      <td>2.34</td>
      <td>Quantitative evaluation of the displacement field in the context of hypomobility: (A) Displacements at the beginning of straining, (B) at the middle of straining and (C) at the end of straining. The white contours are the initial ones.</td>
      <td>111</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Computational flow of the method: The biomechanical model guided by the 3D image-points recovered from intra-abdominal image pair permits to propagate partial tissue deformations to vessels and tumors.</td>
      <td>115</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>3D heterogeneous biomechanical model of the liver with: (a) heterogeneous liver including the vascular network in wire-frame, in (b) beams generated along the vessels.</td>
      <td>118</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Initial model alignment (simulated data): The mesh in purple represents the 3D model computed from CT scans in a different reference, the mesh in blue represents the MLS surface reconstructed from laparoscopic image pair. It represents around 30% of the whole liver that is represented in green.</td>
      <td>119</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>A SURF feature being tracked over time (600 frames) using LK optical flow. (top) tracking performed in both a forward and a backward time direction. Graph shows X and Y coordinates as a function of time. (bottom) 50 × 50 pixel window centered around the feature location at every 100 frames.</td>
      <td>120</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>3.5</td>
      <td>Clustering phase: (left) The two views of 3D point cloud estimated from stereo matching. The Blue sphere represents the rejected features. (right) Image is a representation of the clusters. Red sphere represents the Control Points and the Green sphere the essential features (neighbors).</td>
      <td>121</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>A sequence of images showing a cube being deformed. (Top) augmented reality images where the mesh in wire-frame is superimposed on the video stream, (Middle) the deformed mesh, (Bottom) the volumetric mesh composed of tetrahedra.</td>
      <td>121</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Initial model alignment (simulated data): The mesh in purple represents the 3D model computed from CT scans in a different reference, the mesh in blue represents the MLS surface reconstructed from laparoscopic image pair. It represents around 30% of the whole liver that is represented in green.</td>
      <td>122</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Boundary conditions: The liver model (in green) is constrained by fixed vertices (in purple) that represent the falciform ligament and the underlying fat.</td>
      <td>123</td>
    </tr>
    <tr>
      <td>3.9</td>
      <td>A sequence of images showing the superimposition of the real-time biomechanical model onto the human liver undergoing deformation due to surgical instrument interaction. On the left virtual camera with the original laparoscopic angle and right a different angle.</td>
      <td>124</td>
    </tr>
    <tr>
      <td>3.10</td>
      <td>Volume loss of parotids: (Left) segmentations of the parotids at weeks 1 (red) and 6 (blue). It is worth noticing the volume loss of the parotid as well as the motion of the center of mass. These two parameters have been used to characterize the deformation of parotids in [Bar+04]. (Right) Due to weight loss, parotids may intersect the target volume (in yellow).</td>
      <td>128</td>
    </tr>
    <tr>
      <td>3.11</td>
      <td>Validation: (left) similarity between the initial segmentation and ground truth geometry in blue curve illustrates the deformation of the parotids (DICE decreasing), the red curve exhibits the good similarities between our semi-automatic registration and the ground truth geometry; (middle) planning adaptation using our registration vs no planning adaptation (right). The measured radiation is much lower when the planning is adapted.</td>
      <td>129</td>
    </tr>
    <tr>
      <td>3.12</td>
      <td>Catheter robot from Kings College London with a polymer structure and 4 tendons that enables the robot tip to bend.</td>
      <td>130</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>3.13</td>
      <td>Two canonical cases described where the contacts highly influence the kinematics of the catheter robot. Arrows are the motion directions of the tip for imposed motion at the base: insertion, right cable movement, left cable movement. (a) artificial singularity and (b) inverted mapping result in the deterioration and even the positive-feedback loop.</td>
      <td>133</td>
    </tr>
    <tr>
      <td>3.14</td>
      <td>Framework of decoupled closed-loop bending control for catheter robot through contacts. For the implementation, the FE model of catheter robot and its environment is simulated in real-time. The insertion and bending are decoupled by the projection block which is introduced in Section IV.A. Two independent controllers are employed for the control of catheter robot and its FE model respectively. δ<sub>d</sub> and δ<sub>be</sub> are respectively the desired position and the bending control inputs for the catheter robot, while and δ<sub>FEM</sub> are the corresponding variables for the FE model. δ<sub>in</sub> insertion input for both the catheter robot and the FE model.</td>
      <td>134</td>
    </tr>
    <tr>
      <td>3.15</td>
      <td>Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</td>
      <td>136</td>
    </tr>
    <tr>
      <td>3.16</td>
      <td>Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</td>
      <td>137</td>
    </tr>
    <tr>
      <td>3.17</td>
      <td>Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</td>
      <td>138</td>
    </tr>
    <tr>
      <td>3.18</td>
      <td>Experimental setup for the motion control of a catheter robot through contacts. (a) NDI electromagnetic tracking system: Aurora; (b) a rigid tube for the experiment. The white points are the feature points which are used for the calibration; (c) the cable-driven continuum catheter robot which has four step motors to actuate the bending and one step motor to control the insertion.</td>
      <td>138</td>
    </tr>
  </tbody>
</table>
List of Tables

<table><tbody><tr><td>1.1</td><td>Error measurement. The average and the standard deviation (SD) of the relative energy norm error are given in percentage in the first column. The second column contains the mean displacement of the coil for the three different steps.</td><td>35</td></tr><tr><td>1.2</td><td>CSER computation for two patients. REF corresponds to the CSER of the coil inserted by the clinician, OK represents the CSER of the simulation of the same coil, NOK represent a simulated helical coil. The relative errors to the REF are provided in parentheses.</td><td>51</td></tr><tr><td>2.1</td><td>Evaluation of RBT compared to MHT. SR and SRL are given in percents, TL and CL are given in mm. The mean value was computed over the 10 patients.</td><td>78</td></tr><tr><td>2.2</td><td>Evaluation of LIM accuracy and compacity on 4 synthetic shapes (capsule C, open capsule oC, arc A and bifurcation B) perturbed by a Gaussian noise. d is the 90th percentile of the Euclidean distance of noise free points in P<sub>t</sub> to the output BM. d<sub>a</sub> presents the same statistics using the approximate distance proposed by [Tau91]. N<sub>b</sub> is the number of blobs in the output BM. Minimum (min), median (med) and maximum (max) values are reported over 300 tests.</td><td>82</td></tr><tr><td>2.3</td><td>Distribution in % of the BMs according to the error of fit d<sub>a</sub> in 4 classes (left column): C<sub>1</sub> corresponds to BMs with d<sub>a</sub> &lt; 0.5 v<sub>s</sub>, C<sub>2</sub> to 0.5 v<sub>s</sub> ≤ d<sub>a</sub> &lt; v<sub>s</sub>, C<sub>3</sub> to v<sub>s</sub> ≤ d<sub>a</sub> &lt; 2 v<sub>s</sub> and C<sub>4</sub> to 2 v<sub>s</sub> ≤ d<sub>bm</sub>. v<sub>s</sub> is the voxel size. 9 algorithm configurations were investigated depending on parameters t<sub>g</sub> (targeted precision of fit) and N<sub>s</sub> (maximum model complexity). Computations implied 87564 BMs.</td><td>84</td></tr><tr><td>2.4</td><td>Ratio between the number of triangles #Δ in a triangulated surface and the number of blobs N<sub>b</sub> in the local implicit model for 3 values of t<sub>g</sub> (N<sub>s</sub> = 100). Minimum (min), median (med) and maximum (max) values are reported for 10 patients.</td><td>85</td></tr></tbody></table>
<table><tr><td>2.5</td><td>Computation timings for collision detection on various deployment scenarios, using state of the art implementations. For a comparable amount of primitives, LIM outperforms triangular meshes for collision detection and makes it easier to achieve interactive simulation.</td><td>88</td></tr><tr><td>2.6</td><td>Summary table for results of the error measures between manual and automatic segmentation – the three sub-tables present results of each organ in different pathological cases.</td><td>99</td></tr></table>