Registration by interactive inverse simulation:
application for adaptive radiotherapy

Eulalie Coevoet, Nick Reynaert, Eric Lartigau, Luis Schiappacasse, Jérémie

Dequidt, Christian Duriez

► To cite this version:

Eulalie Coevoet, Nick Reynaert, Eric Lartigau, Luis Schiappacasse, Jérémie Dequidt, et al. Registration by interactive inverse simulation: application for adaptive radiotherapy. International Journal of Computer Assisted Radiology and Surgery, 2015, 10 (8), pp.1193-1200. 10.1007/s11548-015-1175-4. hal-01261302

HAL Id: hal-01261302

https://hal.science/hal-01261302

Submitted on 25 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Registration by interactive inverse simulation: application for adaptive radiotherapy

Eulalie Coevoet<fnref n="1" />,<fnref n="2" />, Nick Reynaert<fnref n="1" />, Eric Lartigau<fnref n="1" />, Luis Schiappacasse<fnref n="1" />, Jérémie Dequidt<fnref n="2" />, Christian Duriez<fnref n="2" />

Received: date / Accepted: date

## Abstract

**Purpose.** This paper introduces a new methodology for semi-automatic registration of anatomical structure deformations. The contribution is to use an interactive inverse simulation of physics-based deformable model, computed in real-time.

**Methods.** The method relies on non-linear Finite Element Method (FEM) within a constraint-based framework. Given a set of few registered points provided by the user, a real-time optimization adapts the boundary conditions and(/or) some parameters of the FEM in order to obtain the adequate geometrical deformations. To dramatically fasten the process, the method relies on a projection of the model in the space of the optimization variables. In this reduced space, a quadratic programming problem is formulated and solved very quickly. The method is validated with numerical examples for retrieving some unknown parameters such as the Young's modulus and some pressures on the boundaries of the model.

**Results.** The approach is employed in the context of radiotherapy of the neck where weight loss during the 7 weeks of the therapy modifies the volume of the anatomical structures and induces large deformations. Indeed sensitive structures such as the parotid glands may cross the target volume due to these deformations which leads to adverse effects for the patient. We thus apply the approach for the registration of the parotid glands during the radiotherapy of the head and neck cancer.

<footnote n="1">¹ Oscar Lambret Hospital, Lille, France</footnote>

<footnote n="2">² INRIA - University of Lille 1, France</footnote>
**Conclusions.** The results show how the method could be used in a clinical routine and be employed in the planning in order to limit the radiations of these glands.

# 1 Introduction

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 established by using a treatment planning system (TPS) [6] [19] 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 [7] [12]. For instance it is reported [8] that xerostomia (loss of saliva) is due to a higher (*than planned*) exposure of the parotid glands while treating throat cancers.

The first motivation of this work 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 [12] [17]. But the challenge remains 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 [9] 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 [7] [16]. *The Finite Element Method (FEM)* has proven to be a good candidate for *intra-patient* registration problems [3] [1].

The second motivation is to improve its current uses which exhibit several downsides. For instance, registration algorithms are often based on applying forces (or displacements) through the analysis of images [9]. This criticism relates to non-rigid registration methods inspired by snakes [20] [18] that mix different energies in a single system. Image-based energies have no physical foundation (for a mechanical standpoint) but, it will create some boundary forces that will be in equilibrium with physics-based deformable models.

Another issue relates to the quality of the image analysis: some boundary motion may not be captured due to a poor analysis of the image leading to an under-constrained system. Finally FEM-based automatic registration can only be predictive when the elastic properties of the registered structures are known. Most of the times, the elastic properties are roughly estimated, which leads to roughly estimated forces inferred from the medical images.
Consequently, the third and main motivation is to provide a registration
tool that can be driven by the physician (for robustness problem) with a
very simple interface, and with an easy and explicit (i.e. no black-box tool)
control on the parameters that have been used for the registration. For our
application, the registration of the parotid glands (or parotids) should be done
using the parameters that are used in clinical studies: an observed volume loss
and motion of the center of gravity [11]. We emphasize that even if only two
parameters are used, the reported deformations are not affine since parts of
the parotids are constrained by the bones of the jaw which leads to bending
deformations in the whole structure that need to be captured. We also insist
that unlike other registration tools, our application do not rely on the image
pixel grey-values. The registration is obtained as a combination of control
points provided by the radiologist which drive an underlying FEM model of
the considered anatomical structures. The first image set is used to perform
the initial segmentation and the following images are only used to provide
visual feedback for the radiologist.

These motivations have led to the development of a new semi-manual reg-
istration method from which we list the following contributions:

- Our method solves a real-time inverse problem on non-linear FEM which has, to our knowledge, not been addressed before. (Existing real-time inverse methods were limited to linear model [4] [5]). This enables to retrieve some missing parameters of the deformations while doing the registration interactively.

- It provides a good control of the registration result by setting manually few registered points. Yet, using non-linear FEM, complex deformations may be captured even if the dimension of the parameter space is reduced. The deformations are constrained by the optimization process and supervised interactively by the physician given the retrieved values of the parameters.

- This new *patient to patient* registration approach is applied in the context of radiotherapy in order to capture the non-rigid motion of the parotids due to weight loss during a therapy.

- A validation study conducted on 7 patients exhibit results whose quality is comparable with manual segmentation / registration while requiring significant less manpower. The decrease of patient exposure to radiations is also highlighted when using our results for adapting the TPS.

The remainder of the paper is organized as follows: the section 2 details the optimization method that uses real-time non-linear FEM with Quadratic Programming (QP) algorithm as well as a numerical validation and the section 3 provides some results related to the clinical application of radiotherapy.

## 2 Real-time optimization method

This section details the formulation that is used for real-time inverse method
on non-linear FEM (static Saint-Venant Kirchhoff model). It describes how it
can be employed to estimate the external loads, pressures, displacements etc. (in the following we use the generic term of boundary conditions) that lead to a given deformation. To estimate the Young's modulus from an observed deformation, we need to know the intensity of the forces applied on the model (like any other inverse methods [10]). If this latter requirement is not reached (only images of the deformation are available), we can estimate a ratio between parameters of different regions (like in [10]) and obtain the displacement (not the forces) of the boundary conditions. This strategy will be used for the registration of the parotids described in the following section since medical images do not provide force measurements.

## 2.1 FEM-based inverse problem

The static FEM used in this paper accounts for non-linear geometrical deformations and integrates over the structure a Hookean constitutive law. In practice, we use the FEM code available in the framework SOFA [14]. During each step i of the simulation, a linearization of the internal forces is computed:

$$f(\mathbf{x}_i) \approx f(\mathbf{x}_{i-1}) + \mathbf{K}(\mathbf{x}_{i-1})d\mathbf{x} \quad (1)$$

where f provides the volumetric internal stiffness forces at a given position x of the nodes (i.e. degrees of freedom of the discretized model), 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_i - x_{i-1}. The lines and columns that correspond to fixed nodes are removed from the system to improve the condition number of the matrix K (at least one fixed node is required for K to be invertible). Static equilibrium (the sum of external and internal force equals to zero) is sought at each step:

$$-\mathbf{K}(\mathbf{x}_{i-1})d\mathbf{x} = \mathbf{p} + \mathbf{f}(\mathbf{x}_{i-1}) + \mathbf{J}^T\mathbf{\lambda} \quad (2)$$

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

- **boundary multipliers λ_b**: we use these constraints to describe the external forces or the external pressures applied on the boundary conditions that creates the deformation. The nodal force distribution and direction is given by J^T which is updated at each step. λ_b is the unknown intensity of the forces or pressures 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(\mathbf{x}, p + dp) \approx f(\mathbf{x}, p) + (\delta f(\mathbf{x}, p)/\delta p)dp$$

In that case, J^T = δf(\mathbf{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 $\lambda_r$**: set interactively by the user to do a local manual registration on a small number of points<fnref n="1" />. Contrary to a lot of existing registration methods, we do not put any force (therefore no energy) to the association of points or on image pixels so $\lambda_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 [13]: **Step I**, a free configuration $x_{\text{free}}$ of the deformable model is found by solving equation (2) 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 forces 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) $$

$\delta_r$ represents a vector between registered and actual positions of points chosen by the user. Thanks to the projection, this vector is directly linked to the Lagrange multipliers used for boundary conditions and parameters. Then a QP problem is set by minimizing the norm of this vector.

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

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

The solution of the QP will provide the best possible values of $\lambda_b$ and $\lambda_p$ to perform the registration without additional energy. If we have a sufficient number of registration points, the matrix of this QP problem is definite, which leads to unique solution. 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$. However, this unique solution could be valid for only one step in the process (as the value of matrix $W$ is computed with a linearization of the internal forces of the FEM model). To obtain a smooth convergence over the

<footnote n="1">¹ We emphasize that these points can be interpolated between the nodes of the mesh using the FEM shape functions. In that case the value of the shape function will be used to fill the rows of matrix $J$ corresponding to the Lagrange multipliers.</footnote>
steps, the values of min, max and ε can be chosen to limit the jump between
steps.

Thanks to the projection, the size of the QP problem is much smaller than
solving the problem in the motion space of equation 2, allowing to solve this
problem in real-time. In practice we use the QP solver available in the CGAL
Library [2].

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

<figure><img src="image_1.png" /><figcaption>Fig. 1: Numerical validation: (a) Mesh composed of 5159 tetrahedra and 1673 points (b) forward simulation by setting pressures in 4 different cavities (c) inverse simulation by registration of 3 points. Same methodology for Young's modulus estimation under known gravity forces: (d) target points (highlighted in red) after setting 3 different Young's moduli (one color by Young's modulus), (e) the resulting deformation once the modulus have been estimated.</figcaption></figure>

## 2.2 Numerical validation

We present here a preliminary validation of the approach using numerical ex-
amples. The experimental protocol is the following: first we create arbitrary
deformation of a deformable object (in our case the *Stanford Bunny* which is
courtesy of the Stanford Computer Graphics Laboratory) by modifying bound-
ary conditions or model parameters, second three relevant points of the model
are chosen and their position is stored when equilibrium is reached and third
the simulation is restarted without the deformation and the three selected
points are registered using the inverse simulation approach (leading to a per-
fect registration since we add a perfect correspondence between the points in
the undeformed or deformed states). We then compare the difference between
the actual values of boundary conditions or parameters that have been used
and the ones estimated through the inverse method.

A first experiment was conducted by applying different pressures in cavities
that were placed in the Stanford Bunny and the inverse problem objective was
to estimate the pressures that led to the deformation (Fig. 1(b)). Our approach
<table><caption>Table 1: Performance of our approach to retrieve pressures (either positive or negative) applied in cavities of our deformable model. With a perfect registration, the error is below 1% and even with a manual registration the estimated values are close to the real ones (less than 5% of error).</caption><thead><tr><th>Cavity</th><th>Applied Pressure</th><th>Estimated Pressure<br>(with perfect registration)</th><th>Estimated Pressure<br>(with manual registration)</th></tr></thead><tbody><tr><td>#1</td><td>25 kPa</td><td>25.1 kPa</td><td>24.9 kPa</td></tr><tr><td>#2</td><td>-30 kPa</td><td>-29.9 kPa</td><td>-26.9 kPa</td></tr><tr><td>#3</td><td>-35 kPa</td><td>-35.1 kPa</td><td>-33.6 kPa</td></tr><tr><td>#4</td><td>20 kPa</td><td>20 kPa</td><td>19.8 kPa</td></tr></tbody></table>

<table><caption>Table 2: Performance of our approach to retrieve Young's Moduli applied to different parts of our deformable model. At initialization, the Young's Moduli are set with an arbitrary value and with a perfect registration, our approach estimates the Young's Moduli with less than 3% of error.</caption><thead><tr><th>Bunny</th><th>Real Young's Modulus</th><th>Initial Young's Modulus value</th><th>Estimated Young's Modulus</th></tr></thead><tbody><tr><td>Right Ear</td><td>1 kPa</td><td>10 kPa</td><td>0.928 kPa</td></tr><tr><td>Left Ear</td><td>5 kPa</td><td>10 kPa</td><td>4.852 kPa</td></tr><tr><td>Body</td><td>2 kPa</td><td>10 kPa</td><td>2.027 kPa</td></tr></tbody></table>

estimates the pressures with less than 1% error for a perfect registration and less than 5% if we let the user perform a manual registration of the three points. The values are listed in the table 1. The second experiment follows the same methodology and was based on a heterogeneous material composed of three different stiffness (then three different Young's moduli to estimate). In this experiment, the deformation was induced by the softness of the model and gravity forces (Fig. 1(d) and (e)). Our approach estimates the three Young's moduli with less than 1% of error (see the table 2).

## 3 Application to adaptive radiotherapy

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. The deformations of the structures due to the weight loss of the patient are not taken into account. 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 Fig. 2). Consequently parotids may be irradiated more than initially planned leading to xerostomia (saliva loss) for 20% of patients.

Yet, with a robust registration performed just before the therapy, the planning could be adjusted even if the parotids are generally poorly visible in the images. Due to the lack of robustness of the contour detection, automatic reg-
istration is not considered as a valid option by physician: they would have to verify all the slice, one by one. And, in addition, if the algorithm fails, they would have to redo the segmentation manually. This is a typical case of use of our inverse approach: Instead of only relying on a pure automatic algorithm, 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. This manual control of the registration point and the direct feedback given by the algorithm could improve the physician confidence on the results if the registration is successful. This way, the therapist has a direct supervision of the parameters that have been used for the registration.

In the case of our application problem, [11] 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). The deformations observed on Parotids are quite large (more than 30 % the size of the structure) and can not be captured with linear models. Hence, an average reduction in the volume of parotid of 0.19cm³ per day of treatment is found as well as a displacement of the parotid barycenter of 3.1mm at the end of treatment, in the median direction. We can not use these two parameters with an affine registration: the parotids are anatomically constrained by the mandible, so when the gravity center is moving or when there is a volume loss, it creates a deformation of the parotids that can not been captured with affine transformation. Given the displacements observed in the image (more than 30 % the size of the structure), the deformation can neither be captured using a linear deformable model. These two parameters are introduced in our registration method and detailed in the next paragraph. The parotids are modeled with FEM models of about 650 tetrahedra and 200 nodes.

The points that are in contact with the mandible are mechanically fixed to be consistent with medical observations.

<figure><img src="image_2.png" /><figcaption>Fig. 2: 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 [11]. (Right) Due to weight loss, parotids may intersect the target volume (in yellow).</figcaption></figure>
## 3.1 Boundary conditions retrieval

The first parameter (motion of the center of mass) is related to the volume loss of neighboring structures of the parotid leading to a global displacement. A constraint is built using a lumped (diagonal) mass matrix $M$ on the FEM model. We can compute the total mass $m = \text{trace}(M)$ Then, the position of the center of gravity $g$ is computed as a weighted sum of the position $q$ of each vertex of the mesh $g = \sum_{j=0}^{N}(m_{jj}/(m+N))q$. Where $m_{jj}$ is the mass on node $j$. This linear relation provides the construction of the Jacobian: $g = J_g q$. Practically, the motion of the center of mass is only significant along axis $x$ (towards the center of the neck), so optimization is done only on this direction. In our experiments, we found that the results are very similar with and without directions $y$ and $z$. The projection $g_x = J_{g,x}q$ is used. QP formulation allows to constrain the direction of displacement of $g_x$ thanks to a unilateral constraint, so that the optimization never finds a solution where the parotid moves in the wrong direction (no volume increase of neighboring structures has been observed during the treatment).

A second constraint is built to apply a volume loss to the parotid. The constraint is built so that if the tissue has a homogeneous properties, an homogeneous volume loss is observed. We compute a weighted normal at each point of the surface of the mesh that is proportional to the surface area of its neighboring triangles. Then, we apply a geometric diffusion of these normals inside the volume. So, at each point, we have a vector that provides the direction of the constraint. We put these directions in the Jacobian vector of the constraint $J_v$. We introduce a unilateral condition in the QP so that the constraint can only reduce the volume. We can concatenate Jacobian vectors $J_{g,x}$ and $J_v$ to obtain the Jacobian matrix of the boundary condition $J_b$. Sometimes, the volume loss of the parotids is not homogeneous. In that case, we divide the parotid in 3 regions, along the principal geometrical axis of the structure and non-uniform volume loss is computed in the optimization.

## 3.2 Interactive inverse simulation for planning update

Our application starts with the geometrical models of the parotids that have been segmented during the initial planning and a CT image of the patient, after several weeks of therapy. First, an automatic rigid registration between the meshes and the new image is performed using the position of mandible bone. Then, the physician is asked to pick several points on the surface of the mesh and register them on the image (see Fig. 3). As this registration is done on a 2D slice of the 3D image, each registered point creates a 2D constraint. 3D registration is achieved when the user places points on different slices. The inverse simulation starts when the number of registered directions is superior to the number of unknowns. Practically, a maximum of 5 values are retrieved during the optimization, so 3 registered points (since each register point induces two constraints) are sufficient. Obviously to improve the precision of the
registration, more parameters or boundary condition values can be optimized,
but the user will have to register more points.

3.3 Validation by comparison with ground truth segmentation

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 achiev-
ing 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 be-
tween 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 com-
puted 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 in
Fig. 4 (left) illustrates that the parotids deformation is significant and second
that our method exhibits good similarity compared to manual segmentation
(average DICE between [0.8;0.9]).

3.4 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

<figure><img src="image_3.png" /><figcaption>Fig. 3: Registration of the parotid deformations: (left) user interface that allows to select 2D points to be registered. (middle) in purple points to be registered on the targeted points (blue). (right) parotid deformation after our inverse simulation.</figcaption></figure>
irradiation map of the parotids without considering planning adaptation Fig.4 (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.

## Conclusion

In this paper we presented a new method for patient to patient registration, based on inverse real-time simulation and an interactive manual registration of a few set of points. This method allows for a good control of the physician on the registration results, which is critical in applications such as radiotherapy. The method could have many other applications including parameter estimation for biomechanics of soft tissues or physics based registration. In future work, the results provided by the inverse method will be confronted to real measurements (for the Young's Modulus for instance) and it will be extended to other elastic parameters (Poisson's ratio). For the radiotherapy application, we will extend the approach to the registration of all the structures around the tumor, instead of only considering the parotids. Moreover, we will investigate more quantitatively how much the therapist can compensate the bad quality of images from cone beam CT (routinely acquired before the treatment) by her/his knowledge and interpretation using our method. We will also look to define the confidence in the image pixels in order to use the richer information given by the image to provide additional constraints to drive the parotids deformation.

## Acknowledgements

Authors would like to thank Pierre Jannin for his advice on validation. The project had the financial support of ANR JCJC Simi3 (Ideas), Oscar Lambret Hospital and Inria Lille Nord-Europe research centre.

<figure><img src="image_4.png" /><figcaption>Fig. 4: 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>
## Conflict of interest

Eulalie Coevoet, Nick Reynaert, Eric Lartigau, Luis Schiappacasse, Jérémie Dequidt and Christian Duriez declare that they have no conflict of interest.

## References

1. Schnabel Julia, Tanner Christine, Castellano-Smith Andy, Degenhard Andreas, Leach Martin, Hose, Rodney Hill Derek and Hawkes David: Validation of nonrigid image registration using finite-element methods: application to breast MR images, IEEE Transactions on Medical Imaging 22(2), 238-247 (2003)

2. Fogel, Efi and Teillaud, Monique: The computational geometry algorithms library CGAL, ACM Communications in Computer Algebra, 47(3/4):85-87(2014)

3. Zhu, Yanning and Hall, Timothy J and Jiang, Jingfeng: A finite-element approach for Young's modulus reconstruction, IEEE Transactions on Medical Imaging, 22(7), 890-901 (2003)

4. Becker, Markus and Teschner, Matthias: Robust and Efficient Estimation of Elasticity Parameters using the linear Finite Element Method, SimVis, 15-28 (2007)

5. Eskandari Hani, Salcudean Septimiu, Rohling, Robert and Bell, Ian, Real-time solution of the finite element inverse problem of viscoelasticity, Inverse Problems journal, 27(8) 85-102 (2011).

6. Webb, S.: The physical basis of IMRT and inverse planning. The British Journal of Radiology, 76, 678-689 (2003).

7. Nelms, B. and Tomé, W. and Robinson, G. and Wheller, J.: Variations in the contouring of organs at risk: test case from a patient with oropharyngeal cancer. Intl Journal of Radiation Oncology,Biology,Physics, 82, 368-378 (2012).

8. Dirix, P. and Nuyts, S. and Van den Bogaert W.: Radiation-induced xerostomia in patients with head and neck cancer. Cancer, 107(11):2525-34 (2006).

9. Crum, W.R. and Hartkens, T. and Hill, D.L.G.: Non-rigid image registration: theory and practice. British Institute of Radiology (2014).

10. Lee H-P, Foskey M., Niethammer M., Krajcevski, P., Lin M.: Simulation-based joint estimation of body deformation and elasticity parameters for medical image analysis. IEEE Transactions on Medical Imaging, 31, 2156-2168 (2012).

11. Baker JL et al., Quantification of volumetric and geometric changes occurring during fractionated radiotherapy for head-and-neck cancer using an integrated CT/linear accelerator system, Intl Jnl of Rad. Oncology Biology Physics (2004).

12. Veiga, C. et al, Toward adaptive radiotherapy for head and neck patients: Feasibility study on using CT-to-CBCT deformable registration for dose of the day calculations. Medical physics, 41(3), 031703. (2014)

13. Guebert, C. Duriez, C. and Grisoni, L.: Unified processing of constraints for interactive simulation. VRIPHYS (2008)

14. Faure, François et al: Sofa: A multi-model framework for interactive physical simulation. In the book "Soft Tissue Biomechanical Modeling for Computer Assisted Surgery", p 283-321 (2012)

15. Baraff, David and Witkin, Andrew: Large Steps in Cloth Simulation. SIGGRAPH '98, 43-54, 1998.

16. Indra J. Das and Vadim Moskvin and Peter A. Johnstone: Analysis of Treatment Planning Time Among Systems and Planners for Intensity-Modulated Radiation Therapy. Journal of the American College of Radiology, 6(7), 514-517, 2009.

17. Pierre Castadot and John Aldo Lee and Adriane Parraga and Xavier Geets and Benoît Macq and Vincent Grégoire. Comparison of 12 deformable registration strategies in adaptive radiation therapy for the treatment of head and neck tumors. Radiotherapy and Oncology, 89(1), 1-12, 2008.

18. Nikos Paragios and Mikael Rousson and Visvanathan Ramesh. Non-rigid registration using distance functions. Computer Vision and Image Understanding, 89(1), 142-165, 2003.
19. Michael Goitein and Mark Abrams. Multi-dimensional treatment planning: I. Delin-
eation of anatomy. International Journal of Radiation Oncology*Biology*Physics, 9(6),
777-787, 1983.

20. M. Kass and A. Witkin and D.Terzopoulos, Snake: Active contour model, Volume 1, (1987)