UNIVERSITE PARIS-SUD

ÉCOLE DOCTORALE : Modélisation et Instrumentation en Physique,

Énergies, Géosciences et Environnement
Laboratoire d’Analyse de Données et Intelligence des Systèmes

DISCIPLINE : Physique

THÈSE DE DOCTORAT

soutenue le 13/03/2014

par

Mingchao YANG

Optimisation des plans de traitement en

radiothérapie grâce aux dernières techniques de
calcul de dose rapide

Directeur de thèse : Jean-Marc BORDY Chercheur (CEA LIST, Saclay)

Composition du jury :

Président du jury : Isabelle BERRY Professeur (Hopital de Rangueil, Toulouse)

Rapporteurs : Véronique DEDIEU Chercheur (Centre jean perrin, Clermont-Ferrand)
Jean-michel LETANG Chercheur (Creatis, Lyon)
Examinateurs : Jean-Marc BORDY Chercheur (CEA LIST, Saclay)
Bruno ESPAGNON Professeur (Université Paris-Sud IPN, Orsay)
Francesc SALVAT Professeur (Université de Barcelona, Barcelona)
Membres invités : David MERCIER Ingénieur-Chercheur (CEA LIST, Saclay)
2
Résumé
Cette thèse s'inscrit dans la perspective des traitements de radiothérapie en insistant sur la nécessité de
disposer d’un logiciel de planification de traitement (TPS) rapide et fiable. Le TPS est composé d'un
algorithme de calcul de dose et d’une méthode d’optimisation. L'objectif est de planifier le traitement
afin de délivrer la dose à la tumeur tout en sauvegardant les tissus sains et sensibles environnant.

La planification des traitements consiste à déterminer les paramètres d’irradiation les mieux adaptés au
patient. Dans le cadre de cette thèse, les paramètres d'un traitement par RCMI (Radiothérapie
Conformationnelle avec Modulation d'Intensité) sont la position de la source, les orientations des
faisceaux et, pour chaque faisceau composé de faisceaux élémentaires, la fluence de ces derniers. La
fonction objectif est multicritère en associant des contraintes linéaires.

L’objectif de la thèse est de démontrer la faisabilité d'une méthode d'optimisation du plan de

traitement fondée sur la technique de calcul de dose rapide développée par (Blanpain, 2009). Cette
technique s’appuie sur un fantôme segmenté en mailles homogènes. Le calcul de dose s’effectue en
deux étapes. La première étape concerne les mailles : les projections et pondérations y sont
paramétrées en fonction de critères physiques et géométriques. La seconde étape concerne les voxels:
la dose y est calculée en évaluant les fonctions préalablement associées à leur maille.

Une reformulation de cette technique permet d’aborder le problème d’optimisation par la méthode de
descente de gradient. L’optimisation en continu des paramètres du traitement devient envisageable.
Les résultats obtenus dans le cadre de cette thèse ouvrent de nombreuses perspectives dans le domaine
de l’optimisation des plans de traitement en radiothérapie.

Mots clés : radiothérapie, système de planification de traitement, TPS, calcul de dose, maille
homogène, descente de gradient

Abstract
This thesis deals with the radiotherapy treatments planning issue which need a fast and reliable
treatment planning system (TPS). The TPS is composed of a dose calculation algorithm and an
optimization method. The objective is to design a plan to deliver the dose to the tumor while
preserving the surrounding healthy and sensitive tissues.

The treatment planning aims to determine the best suited radiation parameters for each patient’s
treatment. In this thesis, the parameters of treatment with IMRT (Intensity modulated radiation
therapy) are the beam angle and the beam intensity. The objective function is multicritiria with linear
constraints.

The main objective of this thesis is to demonstrate the feasibility of a treatment planning optimization
method based on a fast dose-calculation technique developed by (Blanpain, 2009). This technique
proposes to compute the dose by segmenting the patient’s phantom into homogeneous meshes. The
dose computation is divided into two steps. The first step impacts the meshes: projections and weights
are set according to physical and geometrical criteria. The second step impacts the voxels: the dose is
computed by evaluating the functions previously associated to their mesh.

A reformulation of this technique makes possible to solve the optimization problem by the gradient
descent algorithm. The main advantage of this method is that the beam angle parameters could be
optimized continuously in 3 dimensions. The obtained results in this thesis offer many opportunities in
the field of radiotherapy treatment planning optimization.

Key words: radiotherapy, treatments planning system, TPS, dose calculation, homogeneous mesh,
gradient descent

3
4
Remerciements
Les travaux présentés dans ce manuscrit ont été réalisés au sein du laboratoire LADIS (Laboratoire
Analyse de Données et l’Intelligence des Systèmes) du CEA LIST. Tant de chemin parcouru depuis
mes premiers pas sur le territoire français qui m’ont amenés à la réalisation et l’aboutissement de cette
thèse de doctorat. Cette route n’a pu en effet être parcourue que grâce à ceux qui m’ont accompagné et
m’accompagnent encore aujourd’hui. Je souhaite, à présent, leur exprimer toute ma gratitude et mes
remerciements.

Je tiens tout d'abord à exprimer ma profonde reconnaissance à Jean-Marc Bordy, mon directeur de
thèse, pour avoir dirigé ce travail de recherche pendant ces quatre années. Je le remercie pour la
confiance qu’il m’a accordée malgré les périodes de doute, de m'avoir guidé, conseillé et le temps
qu’il a consacré à la relecture de ce manuscrit.

Véronique Dedieu et Jean-Michel Létang m’ont fait honneur d’avoir accepté la lourde tâche de
rapporteur. Je les remercie pour l’intérêt qu’ils ont porté à ce travail, de même que pour leur
participation au Jury. La version finale de ce mémoire a bénéficié de leur lecture très attentive, de leurs
remarques ainsi que de leurs suggestions pertinentes et précieuses. Je leur en suis très reconnaissante.

Je tiens à remercier sincèrement l’ensemble des membres du jury pour avoir accepté d’examiner ce
travail, et qui m’ont fait l’honneur d’assister à sa présentation.

Mes remerciements s’adressent ensuite tout naturellement à mon encadrant David Mercier, pour son
encadrement, ses conseils, sa rigueur scientifique et ses connaissances qui ont permis d’accomplir ce
travail, ainsi que pour ses encouragements durant des périodes critiques et difficiles qui m’ont permis
de ne jamais dévier de mon objectif final. Je le remercie tout particulièrement pour le temps qu’il a
consacré aux nombreuses relectures et corrections du manuscrit, pour les fructueuses discussions que
nous avons eues ensemble et enfin, pour ses compétences dans le domaine de l’informatique. Qu’il
trouve dans l’accomplissement de ce travail toute la sympathie et le respect que je lui porte.

Je remercie également Anthony Larue, actuel responsable du LADIS, pour m’avoir accueillie au sein
de son équipe et avoir mis en œuvre tout ce qui lui était possible pour que cette thèse se termine dans
les meilleures conditions.

Je ne pourrais passer sous silence l’aide et la disponibilité de Baptiste Blanpain, ancien thésard
travaillant sur Doséclair. J’ai pu bénéficier de ses conseils techniques sur la méthode de calcul de dose
et de son précieux encouragement.

Je souhaite remercier aussi Cindy Leloirec pour sa gentillesse et son aide sur le calcul de dose en
Penelope.

Tous mes remerciements à mes collègues du LADIS qui par leur bonne humeur, leur disponibilité et
leur aide m’ont non seulement permis de mener cette thèse à bien mais également de profiter d’une
agréable ambiance de travail. Je pense ici en particulier à Michaël Aupetit, Jean-Philippe Poli,
Laurence Boudet, Lorène Allano, Cédric Auliac, Marine Depecker, Krystyna Biletska, Javier Quijano,
Jérôme Gauthier, Pierre Blanchart, Nicolas Heulot, Maxime Maillot, Khrystyna Kyrgyzova, Antoine
Lachaud, Credo Paniah, Néhémy Lim…

Je dédie ce travail et tout ce qu’il représente :

A ma famille du côté de l’océan pacifique qui étant aussi loin m’a soutenue de si près …

A mes parents, leurs présences et leurs encouragements sont pour moi les piliers fondateurs de
ce que je suis et de ce que je fais.

5
 A ma grande mère.

A ma petite famille en France :

A mon cher époux, sa présence et son soutien quotidien occupent une grande place dans
l’accomplissement de ces travaux. Merci de m’avoir tenu la main jusqu’à la dernière ligne de ce
manuscrit. Merci d’avoir toujours été là pour moi.

A ma petite Prisca et mon petit Samuel que j'aime plus que tout, mes adorables enfants, leur
arrivée nous a procuré la plus grande joie...

Merci à tous pour m’avoir conduit à ce jour mémorable.

6
Table de matières
Glossaire ................................................................................................................................................ 10
Introduction ........................................................................................................................................... 11

Première partie – Contexte scientifique

Chapitre I Nature et objet de l’étude .............................................................................................. 16
1 Le cancer et son traitement ............................................................................................................ 16
2 Définition des volumes en radiothérapie ....................................................................................... 16
2.1 Volumes à traiter ou volumes cibles ..................................................................................... 16
2.2 Volumes relatifs à la dose...................................................................................................... 17
2.3 Volumes à protéger ............................................................................................................... 17
3 Les techniques de la radiothérapie ................................................................................................ 18
3.1 La radiothérapie conventionnelle et la radiothérapie conformationnelle .............................. 21
3.2 La radiothérapie conformationnelle avec modulation d’intensité ......................................... 22
3.2.1 Le principe de la RCMI ................................................................................................. 22
3.2.2 La RCMI par faisceaux stationnaires ............................................................................ 24
3.2.3 La RCMI par faisceaux mobiles en rotation (arcthérapie) ............................................ 25
3.3 Cyberknife ............................................................................................................................. 28
4 Le processus de traitement en clinique .......................................................................................... 30
5 La qualité d’un plan de traitement ................................................................................................. 31

Deuxième partie – Calcul de dose

Chapitre I Dosimétrie ..................................................................................................................... 35
1 Transport de l’énergie par les photons et ionisation ...................................................................... 35
2 Le dépôt de l’énergie par les électrons .......................................................................................... 37
3 Terma, Kerma et dose absorbée .................................................................................................... 38
Chapitre II Les méthodes de calcul de dose : un état de l’art .......................................................... 41
1 Méthode Monte Carlo ................................................................................................................... 42
2 Méthodes basées sur la séparation des rayonnements primaires et diffusés ................................. 42
3 Méthodes de convolution / superposition ...................................................................................... 42
3.1 Méthode point Kernel ............................................................................................................ 43
3.2 Méthode Pencil Beam ........................................................................................................... 44
3.3 Méthode Collapsed cone convolution ................................................................................... 45
4 Méthodes de calcul par réseaux de neurones................................................................................. 46
Chapitre III Représentation analytique de la dose - Doséclair .......................................................... 48
1 Le maillage du fantôme ................................................................................................................. 49
2 Les modèles pré-calculés............................................................................................................... 50
3 Découpage du faisceau primaire en SBIMs .................................................................................. 50

7
4 Fonction de projection ................................................................................................................... 52
5 Calcul de la dose par somme pondérée de projections .................................................................. 53
6 L’algorithme de calcul de dose...................................................................................................... 58
Chapitre IV Contribution Doséclair .................................................................................................. 60
1 Les projections .............................................................................................................................. 60
1.1. Les repères de coordonnées ................................................................................................... 60
1.2. Le point origine du faisceau RCMI exprimé en coordonnée sphérique ................................ 61
1.3. Formulation d’une projection ................................................................................................ 62
2 Propagation axiale ......................................................................................................................... 64
2.1 Tenseur de la maille support.................................................................................................. 64
2.2 Simulation et résultat ............................................................................................................. 66
3 Propagation latérale ....................................................................................................................... 68
3.1 Tenseur de la maille voisine .................................................................................................. 68
3.2 Simulation et résultat ............................................................................................................. 70
4 Discussion ..................................................................................................................................... 72

Troisième partie – Optimisation du plan de traitement

Chapitre I Méthodes d’optimisation de traitement : un état de l’art ............................................... 77
1 Les critères d’optimisation ............................................................................................................ 78
2 Formulation mathématique du problème ....................................................................................... 78
3 Optimisation des intensités seules ................................................................................................. 80
3.1 Modèle linéaire ...................................................................................................................... 80
3.2 Modèle non linéaire ............................................................................................................... 80
4 Optimisation des incidences .......................................................................................................... 81
4.1 Méthodes exhaustives............................................................................................................ 82
4.2 Méthodes itératives déterministes ......................................................................................... 82
4.3 Méthodes itératives stochastiques ......................................................................................... 83
Chapitre II Optimisation des plans de traitement en RCMI............................................................. 85
1 Modèle optimisation et formulation mathématique....................................................................... 85
1.1 Objectifs d’optimisation ........................................................................................................ 85
1.2 La modulation d’intensité du faisceau ................................................................................... 86
1.3 Balistique du traitement......................................................................................................... 87
2 Gradient de dose ............................................................................................................................ 87
2.1 Composition du gradient de dose .......................................................................................... 87
2.2 Le gradient du tenseur de projection ..................................................................................... 90
3 Descente de gradient ..................................................................................................................... 90
3.1 Algorithme générale .............................................................................................................. 90
Chapitre III Implémentation logicielle et validation des résultats .................................................... 92
1 Implémentation logicielle .............................................................................................................. 92

8
2 Validation des paramètres individuels........................................................................................... 93
2.1 Conditions de simulation ....................................................................................................... 93
2.2 Validation du positionnement du point d’origine du faisceau ............................................... 94
2.3 Validation de l’orientation du faisceau.................................................................................. 96
2.3.1 Paramètre ................................................................................................................... 97
2.3.2 Paramètre ................................................................................................................. 100
2.3.3 Paramètre ................................................................................................................. 103
3 Fantôme homogène d’eau avec cible de forme concave ............................................................. 105
3.1 Conditions de simulation ..................................................................................................... 105
3.2 Cyberknife - faisceaux élémentaires indépendants ............................................................. 106
3.2.1 Intensité uniforme (49 mins) ....................................................................................... 106
3.2.2 Intensité variables (24.5 mins) .................................................................................... 110
3.3 Faisceaux RCMI classique .................................................................................................. 114
3.3.1 Plan angles fixes avec modulation (3.7 mins) ............................................................. 114
3.3.2 Plan isocentre (23 mins) .............................................................................................. 117
3.3.3 Plan 3D libre (76 mins) ............................................................................................... 118
4 Fantôme hétérogène Penelope Poumon-Os ................................................................................. 119
4.1 Conditions de simulation ..................................................................................................... 119
4.2 Échantillonnage des points de contrôles.............................................................................. 120
4.3 Étude de cas (82 mins) ........................................................................................................ 120
Conclusion et perspectives .................................................................................................................. 124
Annexe ................................................................................................................................................ 126
Bibliographie ....................................................................................................................................... 152
Publications et communications .......................................................................................................... 156

9
Glossaire
- CTV (Clinical Target Volume) : tumeur visible complétée par une extension tumorale
supposée;

- Fantôme : représentation numérique d’un patient;

- Fantôme homogène : fantôme composé d’un unique matériau;

- Fantôme hétérogène : fantôme composé d’au moins deux matériaux différents;

- Fluence primaire : nombre de photons passant par une unité d’aire;

- Dose : énergie déposée, rapportée à la masse;

- GTV (Gross Tumor Volume) : tumeur visible sur l’imagerie;

- HDV : histogramme dose-volume;

- Index delta, index gamma : techniques d’évaluation de la qualité d’une distribution de dose,
par comparaison à une distribution de référence;

- Interface : surface séparant deux régions de compositions différentes ;

- Isodoses : lignes constituées de points qui reçoivent la même dose;

- IMRT : Intensity Modulated Radiation Therapy (RCMI en Français);

- OAR : organe à risque;

- Projection : association d’un point d’un fantôme maillé à un point d’un fantôme homogène,
les deux fantômes étant irradiés par le même faisceau;

- PTV (Planned Target Volume) : CTV complété par une marge de sécurité tenant compte du
mouvement du patient;

- RCMI : Radiothérapie Conformationnelle avec Modulation d’Intensité (IMRT en anglais);

- SBIM : Sub Beam In a Mesh;

- Scaling : dilatation ou compression d’une distribution de dose, ou du parcours d’une particule,

en fonction de la densité électronique du milieu;

- Tenseur : outil mathématique d’un espace à plusieurs dimensions permettant d’effectuer des
changements de repère;

- Terma : Total Energy Released per unit Mass;

- TPS : Treatment Planning System (système de planification des traitements);

- Voxel : volume élémentaire dans une structure tridimensionnelle numérisée (”voxel” est la
contraction de ”volumetric pixel”).

10
Introduction
La radiothérapie est une technique de traitement du cancer dont un des principes est de détruire les
cellules tumorales. La destruction de la tumeur est obtenue en radiothérapie externes au moyen
d’irradiations successives par des faisceaux de rayonnements ionisants, tout en préservant au
maximum les tissus sains. Le mode opératoire nécessite une connaissance morphologique de la tumeur
ainsi que sa localisation précise dans l’organisme du patient. Ces informations sont obtenues à partir
d’images issues d’un scanner ou encore d’images obtenues par résonance magnétique (IRM).
L’objectif thérapeutique est défini par un radiothérapeute en indiquant pour chaque zone à traiter
(tumeur) une dose minimale ou une dose maximale pour les organes sains.

Différentes techniques de traitement externe existent, parmi lesquelles la Radiothérapie

Conformationnelle par Modulation d’Intensité, ou RCMI, est une des plus complexes à mettre en
œuvre. Le principe de la RCMI est d’adapter la forme du faisceau à celle de la tumeur tout en faisant
varier la fluence de rayonnements au sein du faisceau et de fractionner la dose totale en plusieurs fois,
sous différents angles. La fluence s’exprime en termes de nombre de particules par centimètre carré et
correspond à l’intensité du faisceau. La dose absorbée représente la quantité d’énergie absorbée par
unité de masse. La dose absorbée est directement liée à l’impact du traitement en matière de
destruction de cellules tumorales. La dose absorbée dans l’eau est la grandeur de référence utilisée en
radiothérapie.

La notion de TPS (système de planification des traitements) est apparue avec le développement de
l’imagerie et des codes de calcul en radiothérapie. Cet outil informatique est particulièrement utile
dans le cadre de la RCMI pour utiliser avantageusement les nombreuses possibilités des protocoles
d’irradiation. Le TPS aide tout au long de la chaine de traitement : acquisition des données, délinéation
des structures, définition des faisceaux, calcul de dose et contrôle de la concordance calcul/mesure.
Dans le cadre de cette thèse, nous ne nous sommes intéressées qu’à deux éléments de cette chaine : la
définition des faisceaux et le calcul de dose.

Idéalement, le TPS devrait définir l’ensemble des paramètres des faisceaux du traitement qui
comprend notamment l’énergie, le nombre, la géométrie et la modulation d’intensité des faisceaux.
Toute la difficulté consiste à choisir les meilleurs paramètres qui permettront d’atteindre une
distribution de dose souhaitée dans le volume cible tout en épargnant les tissus sains dans un temps
cliniquement acceptable.

Selon la complexité du cas de traitement, l’énergie est le nombre de faisceaux sont d’abord fixés par le
radiothérapeute. Ensuite, le choix des paramètres restants est effectué en deux étapes :

1. L’optimisation de la géométrie qui consiste à définir les paramètres géométriques des faisceaux tels
que positions et orientations. Dans la majorité des cas, des faisceaux équi-espacées sont utilisés. Ils
sont sélectionnés manuellement par l’utilisateur (radiothérapeute ou physicien) qui raisonne
intuitivement à partir de critères géométriques à l’aide d’un outil du TPS, le « Beam’s eye view »
(BEV). Des méthodes automatiques existent aussi dans la littérature, mais elles sont souvent limitées à
la détermination d’angles coplanaires et discrets.

2. L’optimisation des fluences qui consiste à calculer la modulation d’intensité pour chacun des
faisceaux retenu à l’étape précédente.

Cette démarche est séquentielle. L’inconvénient de celle-ci est que, dans les TPS commercialisés,
l’orientation du faisceau est déterminée intuitivement en fonction des critères géométriques avant de
passer à l’étape d’optimisation des intensités. La modulation des faisceaux ne peut pas être prise en
compte simultanément pour déterminer les orientations optimales. Si le plan obtenu après
l’optimisation des intensités n’est pas satisfaisant, on droit remodifier la géométrie qui nécessite,

11
parfois, plusieurs aller-retour. Ce processus dépend fortement de l’expérience de l’utilisateur et est
consommateur en temps. Les études (Gaede, et al., 2004)(Stein, et al., 1997)(Craft, 2007)(Pugachev, et
al., 2002) montrent aussi que les meilleurs plans de traitement ne sont pas forcément ceux dont les
orientations sont intuitives. Donc, avec cette stratégie de sélection, on ne peut pas garantir l’obtention
de la solution optimale. Concernant des méthodes automatiques non intuitives proposées dans la
littérature, bien qu’elles intègrent l’information de l’intensité optimisée, le temps de calcul pour
trouver une solution optimale est de l’ordre de quelques heures pour des cas simples et une dizaine
d’heure pour des cas plus complexes.

En outre, pour des cas complexes, un petit nombre de faisceaux ne permet plus d’avoir une
distribution de dose satisfaisante. Actuellement, pour surmonter ce problème, un grand nombre de
faisceaux équi-espacés sont utilisés pendant le traitement. C’est l’idée principale de la technique
IMAT et la Tomotherapie. Néanmoins, il faut noter qu’un grand nombre de faisceaux peuvent avoir
des conséquences indésirables dues à la dispersion des faibles doses dans des volumes de tissus sains
plus importants, ce qui augmente le risque de développer un second cancer. En réalité, pour certaines
tumeurs, un petit nombre de faisceaux avec les orientations correctement optimisées peuvent donner
de meilleurs résultats qu’un grand nombre de faisceaux équi-espacés (Rowbottom, et al., 2001). Le
choix de l’orientation est déterminant pour des traitements ayant moins de cinq faisceaux.(Stein, et al.,
1997)(Soderstrom, et al., 1995).

Des faisceaux non coplanaires, pour certaines tumeurs (ex. tumeur paraspinale (Pugachev, et al.,
2001)) peuvent aussi améliorer la distribution de dose significativement. Il faut noter que l’utilisation
de faisceaux non coplanaires permet un degré de liberté supplémentaire dans le choix de la géométrie
optimale, elle augmente considérablement les temps de calcul (car l’espace des solutions est agrandi).
En pratique clinique de RCMI, les temps de traitement augmente également car le manipulateur est
obliger de rentrer dans la salle de traitement pour la mise en place de chaque faisceau.

En partant de ces observations, pour des cas simples, un plan optimisé sera composé d’un petit nombre
fixe (3-10) de faisceaux coplanaires avec leurs orientations soigneusement optimisées. Ces faisceaux
peuvent être délivrés par la technique de la RCMI statique. Aujourd’hui, le défi est de savoir comment
optimiser les orientations des faisceaux de manière quantitative et non intuitive dans un temps
acceptable en clinique. Nous allons proposer, dans le cadre de cette thèse, une nouvelle méthode qui
permet d’optimiser les orientations des faisceaux en partant d’une initialisation de géométrie
correctement définie.

En outre, pour des cas complexes, notre méthode n’est pas restreinte aux faisceaux coplanaires. Quand
à des faisceaux non-coplanaires, au lieu de les délivrer en RCMI statique, la réalisation peut aussi faire
recours à la technique IMAT ou la radiothérapie stéréotaxique (ex. Cyberknife) sur lesquelles notre
méthode est aussi utilisable. Pour le Cyberknife, l’ensemble des géométries optimales des faisceaux
peut être enregistré dans le robot avant le traitement.

En effet, dans un TPS, le calcul d’un plan optimal passe souvent par une optimisation itérative, à
chaque étape le plan courant est évalué, puis amélioré. L’étape d’évaluation consiste à calculer la dose
en chaque point par le plan de traitement courant, et à déterminer l’écart entre cette dose et la
prescription du radiothérapeute. Le temps d’exécution d’une planification automatique dépend
principalement de celui du calcul de la dose. Il faut donc trouver le bon équilibre entre précision et
temps de calcul pour avoir un bon plan en toute confiance. Les TPS actuels utilisent les méthodes de
type convolution/superposition qui semble être le meilleur compromis temps-précision disponible :
calculs de l’ordre d’une minute et trente secondes, pour des erreurs pouvant aller jusqu’à 8%. Une
autre méthode de calcul de la distribution des doses dites « Monte-Carlo » est la plus précise mais
reste encore trop couteuse en termes de temps d’exécution pour une utilisation opérationnelle
systématique en milieu médical.

La technique de calcul de dose est cruciale dans le processus d’optimisation. Pour répondre aux
besoins de l’optimisation en respectant la contrainte de temps de calcul, on cherche une méthode de

12
calcul rapide qui permet de connaitre la variation de la dose en fonction de la variation des paramètres
du faisceau (géométrie et intensité) sans passer par plusieurs calculs de dose afin d’établir des
différences. À notre connaissance, les méthodes de calcul actuellement disponibles en clinique ne
peuvent pas prédire cette variation de dose.

À partir de 2006, Baptiste Blanpain a proposé une nouvelle approche analytique (Doséclair) pour le
calcul rapide de dose. La méthode a été évaluée et perfectionnée pour atteindre aujourd’hui un
compromis vitesse/précision satisfaisant. L’approche mélange deux concepts originaux pour la
radiothérapie :

1. La première hypothèse est que la modélisation peut être faite par maille. Un ‘modèle’ est appliqué
en chaque maille, et l’ensemble est cohérent via les contraintes aux conditions limites.

2. Le second concept est que la dose déposée par un faisceau en un point, quelle que soit la géométrie,
peut être calculée à partir de la définition du faisceau (intensité, géométrie) et de modèles de
régression appris pour les différents milieux sous l’hypothèse, à chaque fois, d’un milieu homogène.
La principale difficulté a consisté à définir comment le passage au travers des différentes
hétérogénéités affecte les paramètres de ces modèles élémentaires.

Grâce à sa capacité à fournir une formulation analytique de la dose dans des mailles plutôt qu'un calcul
direct, Doséclair permet aussi d’accéder à une formulation analytique du gradient de dose qui sera
utilisée dans la fonction de coût. Plus précisément, une fonction de coût étant définie pour qualifier la
pertinence d'une planification de traitement vis-à-vis d'objectif de dose minimale ou maximale dans les
différents organes. La formulation analytique de la dose et de son gradient donne un ensemble de
propriétés très intéressantes pour l’optimisation, au-delà des gains directs en temps de calcul et en
précision :

1. On peut calculer la dose en un point sans avoir besoin des points adjacents. On peut donc envisager
d’utiliser peu de points du fantôme dans le processus d’optimisation permettant de réduire le temps de
calcul car c’est la prise en compte de tous les voxels qui prend le plus de temps.

2. Il est parfaitement possible de calculer le gradient de la dose en un point par rapport à l’intensité, la
position et l’orientation du faisceau. Cette capacité à calculer le gradient donne accès à des méthodes
d’optimisation jusque-là impossible à utiliser avec les méthodes d’estimation de dose actuelles.

L’objectif de la thèse est donc d’explorer ces propriétés et de démontrer la faisabilité d’une méthode
d’optimisation par descente du gradient en implémentant ces nouveaux algorithmes afin d’atteindre les
objectifs : optimiser la géométrie et les intensités du faisceau en respectant au mieux la prescription et
les contraintes de dose dans un temps de 20 minutes sur un ordinateur de bureau. Notons que la
technique de maillage est d’autant plus bénéfique que le nombre de maille utilisé est faible. La
méthode présentée ici est donc particulièrement bien adaptée à des traitements ayant des grosses zones
homogènes telles que la tumeur de la prostate et du poumon qui nécessite un faible nombre de mailles.
L’exploitation d’une telle approche pour des cas ORL est déjà plus hypothétique.

Ce mémoire est constitué de trois parties. La première présente le contexte scientifique ; la seconde
décrit la physique du dépôt d’énergie et les méthodes de calcul de dose ; la dernière décrit les calculs
de gradient et les résultats de l’optimisation.

13
 Première partie

Contexte scientifique

14
Chapitre I Nature et objet de l’étude .............................................................................................. 16
1 Le cancer et son traitement ............................................................................................................ 16
2 Définition des volumes en radiothérapie ....................................................................................... 16
2.1 Volumes à traiter ou volumes cibles ..................................................................................... 16
2.2 Volumes relatifs à la dose...................................................................................................... 17
2.3 Volumes à protéger ............................................................................................................... 17
3 Les techniques de la radiothérapie ................................................................................................ 18
3.1 La radiothérapie conventionnelle et la radiothérapie conformationnelle .............................. 21
3.2 La radiothérapie conformationnelle avec modulation d’intensité ......................................... 22
3.2.1 Le principe de la RCMI ................................................................................................. 22
3.2.2 La RCMI par faisceaux stationnaires ............................................................................ 24
3.2.3 La RCMI par faisceaux mobiles en rotation (arcthérapie) ............................................ 25
3.3 Cyberknife ............................................................................................................................. 28
4 Le processus de traitement en clinique .......................................................................................... 30
5 La qualité d’un plan de traitement ................................................................................................. 31

15
 Chapitre I Nature et objet de l’étude

1 Le cancer et son traitement

Environ 8 millions de personnes déclarent un cancer chaque année dans le monde entier, et plus de 5
millions de ces personnes décèdent de leur maladie. Selon l'OMS (Organisation mondiale de la Santé),
on estime que 20 millions de nouveaux cas de cancer seront diagnostiqués par an d'ici 2020. Avec le
vieillissement de la population, les maladies cancéreuses vont concerner un nombre croissant de
personnes et représenter une cause fréquente de décès (Löf, 2000), ce qui en fait un problème de santé
publique grave.

En termes simples, le cancer est la transformation de cellules saines en cellules qui se divisent
rapidement et en permanence, c’est-à-dire qui se reproduisent bien au-delà des besoins normaux du
corps. Alors que la plupart des cellules saines se divisent et croissent jusqu'à ce qu'elles rencontrent un
autre tissu ou organe en respectant ainsi les limites anatomiques, les cellules cancéreuses continuent de
se développer allant jusqu’à envahir les tissus voisins et même les détruire.

Un nombre significatif de cancers peuvent être soignés par la chirurgie, la radiothérapie et la

chimiothérapie surtout s'ils sont détectés suffisamment tôt.

La chirurgie a pour objectif de retirer la tumeur et d'évaluer sa gravité et son étendue.

Les ganglions proches de la tumeur sont immédiatement analysés. Une intervention chirurgicale peut
être combinée avec d'autres traitements comme la chimiothérapie et la radiothérapie.

La chimiothérapie est un traitement qui consiste à utiliser des médicaments (en général des substances
toxiques) contre les cellules cancéreuses. La majorité des substances chimio-thérapeutiques
fonctionnent par arrêt de la mitose (division du noyau des cellules), ciblant ainsi préférentiellement les
cellules se divisant très rapidement. Mais les substances sont injectées dans la circulation sanguine, la
chimiothérapie traite donc le corps entier et par conséquent, toutes les cellules en prolifération rapide
sont attaquées. Par exemple les cellules responsables de la pousse des cheveux sont également
affectées, ce qui explique que la perte des cheveux est un effet secondaire courant.

Les cellules en prolifération rapide, telles que les cellules cancéreuses, sont plus sensibles au
rayonnement que les cellules saines. La radiothérapie, qui est sujet de notre étude, utilise cette
propriété en administrant des faisceaux de rayonnement pour irradier les zones cancéreuses.

2 Définition des volumes en radiothérapie

2.1 Volumes à traiter ou volumes cibles
Les progrès de l’imagerie et des systèmes informatiques ont permis de définir plus clairement les
volumes d’intérêt en radiothérapie. Nous allons détailler ici les définitions des volumes venant du
rapport (ICRU 50, 1993). Ils sont schématisés en Figure 1 .

16
 Figure 1 Les volumes d’intérêt

Le volume tumoral macroscopique (Gross Tumor Volume : GTV)

C’est celui qui est visible sur l’imagerie (scanner, IRM). Il recevra la dose la plus forte.

Le volume cible clinique (Clinical Target Volume : CTV)

Il comprend le GTV, ainsi que des tissus avec une probabilité tumorale forte même si cela est non
visibles à l’imagerie. La définition du CTV reste encore subjective pour beaucoup de localisation et est
fondé sur l’expérience et les connaissances de la maladie (atteintes ganglionnaires occultes, par
exemple). La définition du GTV et du CTV constitue une part essentielle de la prescription.

Le volume cible planifié (Planning Target Volume : PTV)

Il comprend le CTV et une marge de sécurité qui permet de prendre en compte les incertitudes de
positionnement, les mouvements éventuels des organes et du patient.

L’ICRU recommande d’optimiser les paramètres de la chaîne de traitement pour homogénéiser le plus
possible la dose à l’intérieur du PTV. Il est recommandé de réaliser une planification de manière à ce
que la dose au PTV se trouve entre 95% et 107% de la dose prescrite.

2.2 Volumes relatifs à la dose

Le volume traité
Il s’agit du volume entouré d’une surface isodose spécifiée par le radiothérapeute, correspondant à un
niveau de dose minimal permettant d’atteindre le but du traitement. Idéalement, ce volume traité
devrait correspondre au volume prévisionnel (PTV).

Le volume irradié
C’est le volume de tissus recevant une dose considérée comme significative vis-à-vis de la tolérance
des tissus sains. On pourra évaluer, par exemple, le volume de l’isodose correspondant à 80%, 50%,
ou 25% de la dose prescrite.

2.3 Volumes à protéger

Les organes à risque (Organs At Risk, OAR) sont des tissus pour lesquels il est crucial de limiter
l’irradiation afin de limiter les effets secondaires. Une attention toute particulière doit être portée à la
distribution de dose aux OARs, essentiellement en raison de l’importance des gradients observés en
bordure de volume cible. Comme nous le verrons dans la troisième partie du mémoire, les contraintes
de dose aux OARs interviennent souvent comme des pénalités dans la fonction de coût à optimiser
pour définir le plan de traitement.

17
Trois classes d’organes à risque ont été définies selon leur niveau de morbidité :

- Morbidité sévère : les organes, susceptibles, en cas de lésions graves, d’entraîner la perte
fonctionnelle totale. Par exemple lésions à la moelle épinière rendant paraplégique, lésion à
rétine ou nerfs optique rendant aveugle etc.

- Morbidité modérée : les organes dont la lésion conduit à une perte fonctionnelle importante.
On y retrouve les glandes salivaires, le cristallin, les oreilles, etc.

- Morbidité transitoire : les organes dont la lésion conduit à une perte fonctionnelle mineure
voire nulle. Par exemple la peau ou les muqueuses.

L’organisation du tissu est importante pour déterminer cette morbidité :

- Une architecture en série correspond à un organe à morbidité sévère car la fonction dépend de
toutes ses sous-unités fonctionnelles. Il peut être représenté par analogie avec les circuits
électroniques en série. La rupture d’un seul composant entraîne la perte totale de fonction de
l’organe. La surdose en un point de cet organe altère donc la fonction de l’organe entier. On
s’intéresse alors à la dose maximale reçue par ce tissu. C’est le cas d’organes comme la moelle
épinière dont la section entraîne une paraplégie en aval.

- Une architecture en parallèle correspond à un organe à morbidité modérée, voire faible ou

transitoire. L’organe est constitué de plusieurs sous-unités fonctionnelles plus ou moins
indépendantes les unes des autres. Ainsi, la perte de fonction de l’organe suite à une
irradiation nécessite la destruction d’un nombre significatif de sous-unités. Si le volume
détruit par l’irradiation est réduit, on évite un retentissement sur l’organe et surtout la qualité
de vie du patient. Ainsi, un tel organe peut recevoir une dose élevée si une partie du volume
est préservée. On s’intéresse donc à une contrainte de type dose moyenne ou dose-volume,
c’est-à-dire qu’une partie du volume ne doit pas être irradiée au-delà d’une certaine dose. On
peut citer comme organe en parallèle les parotides ou encore la rétine.

Finalement, pour chacun des organes en série ou en parallèle, des relations dose-volume doivent être
respectées. Cette relation peut être représentée par des histogrammes dose-volume (HDV), qui est
expliqué dans la section 5 de ce chapitre.

En pratique, la notion du PTV ont été étendues aux OARs. (ICRU 62, 1999) défini un volume
prévisionnel pour les organes à risque (Previsonal risk volume, PRV). Ce volume correspond au
volume des OARs étendu par une marge prenant en compte les mouvements ou les déformations des
OARs à l’intérieur du corps, ainsi que les conséquences des incertitudes de mise en place du patient
durant le traitement. Les PRV sont préférentiellement utilisés pour les organes en série (Aventin, et al.,
2012).

3 Les techniques de la radiothérapie

Un traitement de radiothérapie est personnalisé en fonction de la situation particulière de chaque
patient, qui bénéficiera donc de soins spécifiques adaptés à son cas personnel.

Selon le type de tumeur et sa localisation, différents modes de radiothérapie sont utilisés en clinique,
illustré en Figure 2.

18
 Figure 2 Trois modes de radiothérapie (Coucke, 2013)

La curiethérapie consiste à introduire des substances radioactives dans le corps en les plaçant dans un
espace creux naturel, dans la tumeur elle-même ou à proximité immédiate. L’irradiation est en général
isotrope. Il faut donc bien étudier la répartition des sources pour que la tumeur soit bien détruite tout
en minimisant l’impact sur les tissus sains, ce qui est le principe même de l’optimisation en
curiethérapie.

La radiothérapie métabolique est surtout utilisée dans certaines formes de cancers de la thyroïde.
Dans ce cas, la substance radioactive est administrée par voie orale ou intraveineuse et va se fixer
préférentiellement sur les cellules cancéreuses.

La radiothérapie externe est la plus fréquente. La source de rayonnement est à l'extérieur du patient.
Elle consiste à administrer les rayons à travers la peau et les tissus pour irradier toute la région touchée
par la tumeur ainsi qu’éventuellement les ganglions lymphatiques les plus proches.

La radiothérapie externe utilise un flux de photons suffisamment énergétiques pour interagir avec les
atomes et éjecter les électrons de leur orbite. Ces derniers « voyagent » dans le tissu et déposent leur
énergie. L'énergie déposée par unité de masse de tissu est appelée la dose absorbée. La dose absorbée
est exprimée en gray (Gy). Un gray est égal à un Joule (J) d’énergie absorbée dans un kilogramme (kg)
de matière. Nous détaillerons les principes de la physique dans la deuxième partie du mémoire mais
retenons pour le moment que lors de l’irradiation, les ionisations et l’énergie absorbée provoquent des
dégâts dans les cellules (au niveau de l’ADN entre autres). 2 Gy représente une dose de rayonnement
quotidienne qui est généralement tolérée par les cellules saines. On exploite cette particularité pour la
radiothérapie externe par fractionnement. Par exemple, une dose totale de 60 Gy peut être délivrée par
fractions de 2 Gy en 30 jours de traitement. Pour chaque traitement, la dose prescrite et son
fractionnement dépendent donc de la localisation et de la nature de la maladie.

Dans la suite de ce mémoire, nous allons nous intéresser uniquement à la radiothérapie externe. Le
tableau ci-dessous résume les principales caractéristiques des différentes techniques de la radiothérapie
externe que nous allons présenter en détail dans les sections 3.1, 3.2 et 3.3.

19
 Radiotérapie RCMI RCMI IMAT Tomotherapie Cyberknife
conventionnelle S&S SW
Energie du 4 — 25MV 4— 4— 6MV 6MV 6MV
faisceau 25MV 25MV

Utilisation d’un × × × × Collimation

collimateur fixe
multilame

Déplacement × × ×
continue des
lames pendant
l’irradiation

Irradiation selon × × × × Très

des angles nombreux
discrets

Déplacement × ×
continue du
bras pendant
l’irradiation

Conformation × × × ×
du volume
d’irradiation à
la tumeur

Modulation de × × × ×
l’intensité

Variation du ×
débit de dose du
LINAC

Déplacement du ×
patient (de la
table) pendant
l’irradiation

20
3.1 La radiothérapie conventionnelle et la radiothérapie
conformationnelle

Figure 3 Exemple d’isodose en radiothérapie conventionnelle (Holder, et al., 2004)

En radiothérapie conventionnelle, en règle générale, les faisceaux sont assez larges pour irradier toute
la cible à partir des angles d’irradiation. La planification de traitement se faisait « à la main » d’où
l’utilisation de faisceaux présentant des profils plats et l’utilisation de 2 à 4 faisceaux positionnés sur
les 4 points cardinaux afin de faciliter le calcul déjà complexe. montre un exemple de traitement d’un
cancer du foie avec deux faisceaux parallèles opposés. L'intersection des deux faisceaux crée une zone
de forte dose près d’une forme rectangulaire qui englobe presque tout le volume irradié du patient.
Malheureusement, cette zone contient une structure critique – la moelle épinière. Pour que ce
traitement soit viable, la dose prescrite par le médecin devra être maintenue en dessous de la dose
tolérée par la moelle épinière. Mais la dose peut alors être non conforme dans la tumeur. Cette
approche classique, qui utilise de 2 à 4 faisceaux pour traiter une tumeur, a été la pierre angulaire de la
radiothérapie pendant des années.

Figure 4 Exemple d’isodose en radiothérapie conformationnelle (Holder, et al., 2004)

Avec l’avènement de l’ordinateur et le développement du collimateur multilâmes (MLC), la

radiothérapie conformationnelle (3D-CRT) a fait son apparition à la fin des années 1990. Les scanners
permettent des reconstructions en 3D du corps et de tous les organes. Les logiciels de Beam Eye View
(BEV), vue depuis la cible, permettent de réaliser de façon virtuelle des plans de traitement en 3D qui
contourent plus précisément la tumeur en épargnant les tissus sains. Ceci permet de délivrer une
distribution de dose ayant un très haut degré de conformité avec la forme de la tumeur.

Ces distributions de dose sont représentées par ce que l’on appelle les courbes d'isodoses. Ces
dernières illustrent les iso-niveaux de la dose absorbée. Le niveau de l’isodose est définie par un
pourcentage de la dose cible prescrite. La Figure 4 représente une distribution de dose avec différents
niveaux. La région à forte dose est représentée par la ligne 60 Gy (ligne noire), qui suit la forme de la

21
tumeur. La courbe extérieure est l'isodose à 20%, ce qui signifie que le tissu à l'intérieur de cette
courbe reçoit jusqu’à 20% de la dose prescrite. En délivrant la dose la plus élevée selon la forme
exacte de la tumeur, les tissus sains et critiques à proximité sont épargnés.

Techniquement, pour réaliser cette conformation, les collimateurs multilames sont utilisés pour
remplacer les champs carrés. L’exemple le plus démonstratif est le traitement du cancer de la prostate.
Cet organe de 40 cm3 a la forme d’une pyramide. Avant la fin des années 1990, il est irradié par des
champs carrés de 8 cm de côté délivrant une dose de 60Gy pour ne pas léser le rectum. Ces champs
carrés irradient un volume proche de 500 cm3 (demi-litre) soit dix fois plus environ de tissu sain que
de tissu tumoral. Avec la nouvelle technique de la 3D-CRT, la prostate est irradiée avec des faisceaux
reproduisant la forme de la pyramide et non plus en forme de cube. Le rectum est protégé ainsi que la
vessie, sachant que la dose dans le volume cible a été augmenté à 70-76 Gy. Des essais randomisés ont
montré que cette augmentation de dose s’accompagne d’une amélioration du contrôle local et de la
survie sans augmentation des complications rectales, vésicales voir sexuelles (Pollack, et al., 2000).
Une forte dose dans un petit volume reste un modèle de préférence, sous réserve d’une exactitude du
ciblage balistique. En 2000, la 3D-CRT devient la technique de routine pour la grande majorité des
irradiations notamment à visée curative (tumeurs cérébrales, ORL, poumon, prostate…).

Malgré ce progrès thérapeutique apporté par la 3D-CRT, il faut noter que cette technique est moins
précise dans le cas où la tumeur présente des formes concaves avec en plus des organes à risque dans
ces concavités, ce qui représente 30% des cas. D’où le développement de la Radiothérapie
Conformationnelle par Modulation d’Intensité (RCMI) que nous allons présenter dans la prochaine
section.

3.2 La radiothérapie conformationnelle avec modulation d’intensité

Le développement de la Radiothérapie Conformationnelle par Modulation d’Intensité (RCMI) (IMRT
Collaborative Working Group, 2001) (ANAES, 2003) représente une étape importante dans la
réalisation de la radiothérapie conformationnelle. Elle permet de contourer précisément la tumeur et de
moduler le faisceau pour présenter des isodoses concaves permettant le traitement de la tumeur en
épargnant les tissus sains.

3.2.1 Le principe de la RCMI

Figure 5 Représentation synthétique des travaux de (Brahme, et al., 1982)

22
(Brahme, et al., 1982) a abordé le problème des tumeurs concaves par une illustration symétrique et
synthétique, montré en Figure 5 : une tumeur ronde entourant un organe à risque. En clinique, cette
situation est souvent rencontrée dans les traitements des tumeurs de la tête et cou, dans lesquels les
tumeurs et les ganglions lymphatiques régionaux se situent autour de la moelle épinière. Il a été
montré que, pour la radiothérapie rotationnelle, il est possible de calculer un profil d’intensité
hétérogène pour avoir une distribution de dose homogène dans tel volume cible. Notons que dans la
radiothérapie rotationnelle, le faisceau tourne autour du fantôme. Cette rotation est relative, dans
l’illustration de la Figure 5, elle est modélisée par une rotation de fantôme. Dans ce modèle idéal,
l’atténuation du faisceau n’est pas prise en compte.

On peut constater que le point rentre dans le champs plus tard et sorte du champs plus tôt que le
point . Si le faisceau est de fluence uniforme (Figure 5 (c)), la dose déposé au point est
inférieur à la dose déposée au point . Avec un profil d’intensité hétérogène (Figure 5 (d)),
une égalisé de la dose déposée entre les deux point est obtenue. L’idée de Brahme est démontrée par
(Mackie, 1993). Sur la Figure 6, une distribution de dose homogène dans la cible et une dose nulle au
centre est obtenue à l’aide de 51 faisceaux correctement modulé. L’utilisation des faisceaux avec
modulation d’intensité révèle un potentiel énorme pour le traitement de volumes concaves ou proches
d’organe à risque à éviter. Nous allons détailler les méthodes pour déterminer les intensités optimales
des faisceaux dans la partie III du mémoire.

Figure 6 Une distribution de dose homogène dans la tumeur de forme disque obtenue par des faisceaux d’intensité
modulé (Mackie, 1993)

Maintenant supposons que l’on dispose des intensités optimales, afin de pratiquer cette technique en
clinique, Brahme a proposé d’utiliser des compensateurs métalliques pour moduler les intensités.
Cependant, c’est une technique très lourde car les compensateurs doivent être fabriqués pour chaque
traitement spécifique et ensuite changés manuellement entre chaque faisceau.

23
 Figure 7 Un collimateur mutilâmes de la RCMI classique (Holder, 2005)

La première utilisation d’un collimateur multilames pour la modulation du faisceau est introduite par
la société NOMOS. C’est un collimateur binaire utilisé dans la tomothérapie, nous allons le détailler
dans la section 3.2.3.2.

Aujourd’hui, la plupart des traitements en RCMI sont réalisés à l’aide des collimateurs multilames
montrés en Figure 7. Il s’agit d’un élément mécanique fixé à la tête de l’accélérateur linéaire (Varian
Medical Systems) ou intégré (Elekta, Siemens) qui permet, à l’aide de lames mobiles effectuant des
mouvement continus et asynchrones, de définir précisément la forme du champ et d’atténuer les
rayonnements en dehors du champ.

Grâce à ces collimateurs, la dose est délivrée en deux modes distincts. Mode stationnaire : le bras est
fixé pendant l’irradiation, détaillé dans la section 3.2.2. Mode dynamique : le bras se déplace en
continue pendant l’irradiation, détaillé dans la section 3.2.3.

3.2.2 La RCMI par faisceaux stationnaires

Pour les techniques de RCMI par faisceaux stationnaires, le bras de l’accélérateur reste fixe pendant
l’irradiation. En générale de cinq à douze faisceaux sont utilisés. Ces faisceaux provenant des angles
différents convergent vers un même point. À chaque emplacement fixe, le patient est irradié par des
faisceaux de formes spécifiques (segment) définies par le collimateur, illustré par la Figure 7. La
modulation d’intensité est réalisée par la superposition de segments. Chaque segment possède sa
propre pondération en dose : un certain nombre d’unités moniteurs (UM) est associé à chaque
segment. Selon le mouvement des larmes pendant l’irradiation, deux modes de délivrance des
segments existent : le mode statique nommé Step and shoot (S&S) et le mode dynamique nommé
Sliding window (SW).

3.2.2.1 Mode Step and shoot

En RCMI type S&S, l’irradiation se fait de manière séquentielle. Entre deux irradiations successives,
une nouvelle forme de champ est définie par le déplacement des lames du collimateur. Le champ
d’irradiation est décomposé en plusieurs segments dont le nombre est en fonction des contraintes de
modulation du faisceau. Selon la localisation tumorale à traiter et les contraintes imposées, le nombre
de segments varie entre 2 et 20. En Figure 8, un champ d’irradiation hétérogène (à droite) est
décomposé en 5 segments de formes diverses. Ce mode de délivrance n’est pas limité par la vitesse
des lames, mais il subsiste un temps mort entre les segments qui diminue l'efficacité d'exécution du
traitement. La résolution spatiale de la modulation est alors limitée par la quantité de segments par
incidence et par le nombre d'unités moniteur minimal par segment. Il importe d'avoir un algorithme de
segmentation qui limite la quantité de segments ainsi que le nombre d'unités moniteur total.

24
 Figure 8 Exemple de faisceau possédant 5 segments en technique de RCMI type Step and shoot (Bey, 2010)

3.2.2.2 Mode Sliding window

En RCMI type SW développée par (Convery, et al., 1992) , les lames se déplacent en continu au cours
de l’irradiation, illustré par Figure 9 . Dans ce cas, le faisceau modulé en intensité résulte du balayage
du champ par une fente glissante de forme variable définie par le collimateur mutilâmes. C’est la
variation de vitesse de déplacement des lames qui va créer les différences de fluence pour générer le
faisceau modulé. La résolution spatiale de la modulation est limitée par la vitesse de déplacement des
lames.

Figure 9 Exemple de faisceau modulé délivré en mode Sliding window (Lafond, 2013)

Pour la RCMI par faisceaux stationnaires, le nombre de faisceaux utilisés varie en fonction des
contraintes spécifiques liées à chaque traitement. Dans le cas du traitement de la prostate par RCMI,
cinq faisceaux modulés quasiment équi-répartis permettent une bonne conformation de la dose au
volume cible et une bonne préservation des organes à risque. L’augmentation du nombre de faisceaux
n’apporte pas d’amélioration majeure. Cependant, dans le cas du traitement des tumeurs ORL, pour
obtenir une bonne homogénéité de la dose au volume cible et une bonne protection des OARs, le
nombre de faisceaux nécessaire s’élèvera de sept à neuf. Le nombre de faisceaux est limité dans ce
type de RCMI stationnaire à cause de la contrainte du temps de traitement en clinique. Cela signifie
que le processus d'optimisation doit choisir un certain nombre d’angles optimaux à partir de
l’ensemble des angles possibles. La plupart des logiciels de planification ne procèdent pas à la
sélection automatique des faisceaux, l’utilisateur doit les définir au préalable. En revanche, un tel
logiciel fournit des outils de visualisation qui permettent à l'utilisateur de sélectionner intelligemment
les angles. Cette technique permet de réduire la complexité du problème à un niveau gérable, mais il
ne garantit pas une solution optimale globale.

3.2.3 La RCMI par faisceaux mobiles en rotation (arcthérapie)

3.2.3.1 IMAT

L’arcthérapie avec modulation d’intensité, appelée Intensity Modulated Arc Therapy (IMAT) a été
introduite par (Yu, 1995). Elle est une technique en plein essor, dont l’objectif principal est de
diminuer le temps de traitement et le nombre d’unité moniteur délivrées par rapport à la RCMI
stationnaire tout en maintenant au moins la même qualité de distribution de dose.

C’est une approche dynamique qui fait varier à la fois la position des lames et l'angle de bras de
l'appareil pendant l'irradiation. Pour les cas complexes, le temps de traitement pourra être réduit
considérablement car l’irradiation est continue pendant la séance de traitement.

25
En outre, IMAT présente la flexibilité de conformer du mieux possible le dépôt de dose à la
localisation de la tumeur à l'aide d’un grand nombre d’angles de faisceaux disponibles. (Burman, et
al., 1997) a montré que chaque augmentation du nombre d'angles de faisceaux conduit à une dose plus
homogène dans la tumeur et une diminution de dose aux OAR.

Un arc est défini par deux positions extrêmes entre lesquelles le faisceau est présent tout au long de
l’irradiation. L’arc complet est approximé par une série de champs fixes uniformément espacés que
l’on appelle les points de contrôle (CP). En pratique, chaque CP est espacé de 5 à 10 degré (Yu, 1995)
(Shepard, et al., 2007) (Earl, et al., 2003) (Gladwish, et al., 2007) (Yu, et al., 2002) (Wong, et al.,
2002). Pour chaque CP, une intensité est optimisée à l’aide de la même technique utilisée en RCMI
stationnaire. Ensuite, ces intensités sont converties en plusieurs segments. IMAT doit tenir compte des
restrictions de mouvement des lames du collimateur linéaire lorsque le bras se déplace d'un angle à
l'autre pour réaliser ces segments. En générale, la position des lames du collimateur linéaire est
modifiée de 2 à 4 cm entre deux angles adjacents pour passer d’un segment à l’autre (Yu, 1995)
(Shepard, et al., 2007) (Earl, et al., 2003) (Gladwish, et al., 2007) (Yu, et al., 2002) (Wong, et al.,
2002). Parfois les segments optimisés doivent être modifiés afin de délivrer un segment réalisable.

Pour surmonter cette limitation, l’utilisation des arcs superposés est proposée dans la littérature (Yu,
1995) (Shepard, et al., 2007) (Earl, et al., 2003). Cependant, ces arcs superposés allongent le temps de
traitement.

En outre, l'optimisation s'effectue à l'aide d'un échantillonnage des angles fixes du bras très espacé.
Ceci est susceptible de dégrader la précision du plan, ce qui entraîne potentiellement une erreur
dosimétrique (Shepard, et al., 2007).

En revanche, un échantillonnage trop proche réduit l'amplitude du mouvement des lames entre les
angles adjacents, ce qui empêche algorithme d'optimisation d'exploiter pleinement les avantages d'un
grand nombre de faisceaux.

En 2008, (Otto, 2008) a proposé une nouvelle méthode d’optimisation des plans de traitement pour
IMAT en utilisant plus d’évolution des paramètres de parcours sur un seul arc. Il a alors introduit le
terme VMAT, pour Volumetric Modulated Arc Therapy. En comparant avec IMAT, il est maintenant
possible de faire varier la vitesse de rotation du bras ainsi que le débit de dose de l'appareil pendant
l’irradiation.

Cette méthode utilise un échantillonnage progressif pendant l’optimisation. Les segments sont
optimisés initialement pour un certain nombre d'angles très espacés avec peu de pour la connectivité
entre segments. Pendant l’optimisation, des angles sont ajoutés progressivement jusqu’à l’obtention
d’une certaine fréquence d'échantillonnage prédéfinie. Les contraintes liées à la variation du débit de
dose, le déplacement des lames et la rotation du bras sont pris en compte pendant l’optimisation. La
forme du segment inséré est initialisée par une interpolation linéaire de ses deux segments voisins.
Donc, la connectivité des segments est prise en compte pendant les itérations. L'arc unique optimisé
peut être livré dans les 2 minutes.

Cette technique est connue commercialement sous les appellations VMAT (Elekta, Stockholm,
Suède), RapidArc (Varian Medical Systems, Palo Alto, CA) ou Smart Arc (Philips Healthcare,
Andover, MA).

3.2.3.2 Tomothérapie

La tomothérapie proposé par (Mackie, 1993), se présente sous forme d’un anneau qui tourne autour de
la table de traitement sur laquelle est installé le patient, illustrée en Figure 10. L’anneau est composé
d’un accélérateur linéaire avec un collimateur binaire. Ce collimateur est composé de 64 lames

26
contrôlées par des vannes d'air individuelles, illustré par Figure 11, qui s’ouvrent pour laisser passer le
rayonnement. Avec une vitesse d’ouverture et de fermeture de 17 millisecondes, ces lames mobiles
permettent de cibler de façon très précise la tumeur à atteindre et ainsi de moduler les faisceaux en
fonction du traitement du patient. Un exemple synthétique est illustré en Figure 12, les trois faisceaux
sont modulés de tel sorte que l’OAR est protégé. La table de traitement peut également se déplacer
longitudinalement, ce qui permet d’irradier des très petits champs mais aussi des champs pouvant aller
jusqu’à 1m80 sans s’arrêter, contrairement aux accélérateur classiques qui sont limités à 40cm de
longueur au maximum.

Figure 10 Un traitement en tomothérapie (Shepard, et al., 1999)

Figure 11 Un collimateur mutilâmes de la tomothérapie (Van Dyk, et al., 2002)

Figure 12 Faisceau modulé de tomothérapie (Langen, 2006)

Chaque rotation complète de l'accélérateur linéaire peut être approximée par une série de pas discrets.
Le temps d’ouverture de chaque lame est variable, ce qui affecte l’intensité du faisceau correspondant.
La Figure 13 montre une simulation d’ouverture des lames en deux angles différents. En optimisant
ces intensités des faisceaux, il est possible de générer une distribution de dose proche de la forme de la
tumeur.

27
 Figure 13 Simulation des ouvertures d’un collimateur de la tomothérapie (TomoTherapy, 2011)

3.3 Cyberknife
Le système CyberKnife, développé depuis 1990, est un système de radiothérapie en conditions
stéréotaxiques commercialisé par la société Accuray, Inc., Sunnyvale, CA. Il permet de traiter des
tumeurs de manière non intrusive. L'un des atouts principaux du CyberKnife réside en sa flexibilité; en
effet ce dernier permet d'accéder et de suivre en temps réel des tumeurs présentes sur n'importe quelle
partie du corps, et ce avec une précision de l'ordre du millimètre.

L’irradiation par Cyberknife utilise un accélérateur linéaire miniaturisé monté sur un bras robotisé
piloté par ordinateur, qui permet d’orienter les faisceaux avec six degrés de liberté (antérieur-
postérieur, supérieur-inférieur, gauche-droite, rotation dans les 3 axes (roulis, tangage, lacet)),
explicité par six axes de mobilité illustrés en Figure 14 :

Figure 14 Cyberknife de la société Accuray

Pendant un traitement, la trajectoire du bras est définie par un ensemble de nœuds (possibilité de 120
nœuds), ces derniers sont distribués sur une sphère enveloppant la tumeur. Le centre de la sphère
coïncide avec l’isocentre de la tumeur (isocentre d’imagerie), illustré sur la Figure 15.

28
 Figure 15 Ensemble de nœuds sur la sphère enveloppant la tumeur (Gagnon, 2012)

Pour les tumeurs de forme sphérique, les faisceaux convergent vers cet isocentre pour avoir une
distribution sphérique de haute dose. Sur la Figure 15, un faisceau isocentrique est présenté par le trait
bleu. L’utilisation d’un seul collimateur circulaire de diamètre fixe (même diamètre que la tumeur) est
suffisante. Les douze collimateurs de diamètres différents sont présentés sur la Figure 16.

Figure 16 Collimateurs circulaires fixes

Pour des tumeurs de forme complexe, le traitement est non isocentrique. Plusieurs directions des
faisceaux par nœud sont initialisées aléatoirement. Les intensités du faisceau, les diamètres du
collimateur et les directions du faisceau sont optimisés ultérieurement (Tombropoulos, et al., 1999)
(Schweikard, et al., 2006). Les faisceaux sont présentés en rouge sur la Figure 15. Un exemple de
faisceaux optimisés pour un cancer du poumon est présenté Figure 17.

Figure 17 Cancer de poumon traité par les faisceaux fins de Cyberknife (Colorado Cyberknife, 2014)

Généralement il faut compter entre 100 et 200 faisceaux par traitement (Gagnon, 2012). Le temps de
traitement pour les tumeurs intracrâniens est de 20 – 45minutes et 1 – 2 heures pour les tumeurs
extracrâniennes.

29
Le cyberknife offre une très grande flexibilité de ciblage et d’orientation du faisceau, et permet de
réaliser aisément des traitements en temps réel avec des faisceaux non coplanaires et non
isocentriques. Cependant, le temps de traitement reste à améliorer pour les traitements extracrâniens.
Les tumeurs extracrâniennes possèdent souvent de grandes tailles qui exigent un grand nombre de
faisceaux fins pour assurer le taux de couverture de la tumeur, ce qui rallonge le temps de traitement.
Depuis 2012, un collimateur multilâme (InCise™) est disponible pour le Cyberknife. Ce collimateur
utilise 41 paires de lames en tungstène pour moduler la forme du faisceau non isocentrique et non
coplanaire. Pour un même traitement, le nombre de faisceaux est réduit en utilisant le collimateur
multilâme InCise™, qui implique une réduction de 50% du temps de traitement initial (Accuray,
2013). Ainsi, un logiciel spécifique d’optimisation de trajectoire est mis en place pour augmenter la
rapidité du traitement (Dieterich, et al., 2011).

Nous allons voir dans la troisième partie du mémoire le potentiel de notre méthode appliquée sur cette
technique de radiothérapie.

4 Le processus de traitement en clinique

Un cancer est souvent diagnostiqué suite à l'observation de symptômes liés à la maladie chez un
patient. Le patient est alors désigné pour les études d'imagerie qui peuvent être la tomodensitométrie
(TDM ou CT pour Computed Tomography), l’imagerie par résonance magnétique (IRM) ou la
tomographie par émission de positons (TEP).

Pour la RCMI, le corps médical utilise un logiciel appelé Système de Planification des Traitements
(TPS) qui est un outil informatique servant à assister le médecin et le physicien pour la définition et la
réalisation du traitement.

Le processus de traitement commence par la construction d’un modèle 3D du patient à partir des
multiples coupes issues d’une ou de plusieurs des méthodes d’imageries mentionnées précédemment.
Le TPS offre ensuite des outils de contourage qui permettent de délimiter les tumeurs et les organes à
protéger. Les données anatomiques (contour de tumeur et des organes à protéger) servent à la
définition des volumes d’intérêt tels que les PTVs et OARs. La dose prescrite pour chaque PTVs et les
contraintes de dose pour les PTVs et les OARs sont définies. La démarche classique de la définition
des faisceaux d’irradiation tel que leurs caractéristiques physiques (énergie et type) et géométriques
(nombre, incidences…) est basée sur l’expérience et l’intuition du corps médical.

Une fois les faisceaux sélectionnés, le TPS est utilisé pour optimiser la carte des fluences pour chaque
faisceau. La puissance informatique moderne permet de calculer des cartes de fluences optimales en
20 minutes environ.

Pendant ce processus d’optimisation, le TPS va calculer la distribution de dose déposée par les
faisceaux dans les volumes anatomiques. Ces volumes anatomiques sont généralement découpés en
petits éléments de volume réguliers (voxels). Donc, ce calcul consiste à modéliser la dose déposée
dans chaque voxel.

Si la distribution finale de dose déposée générée par les différents faisceaux sélectionnés n’est pas
acceptable, le processus de sélection des faisceaux et d’optimisation de la carte de fluences se répète
jusqu'à ce qu'un traitement satisfaisant soit trouvé mais les temps de calcul sont alors plus importants.
Nous allons détailler les méthodes de calcul de dose dans la deuxième partie du mémoire et le
processus d’optimisation dans la troisième partie du mémoire.

30
5 La qualité d’un plan de traitement
Afin de vérifier que le plan de traitement proposé par le TPS satisfait au mieux les contraintes
médicales préétablies, des évaluations qualitatives et quantitatives de la distribution de dose sont
nécessaires.

D’un point de vue qualitatif, l’évaluation des courbes isodoses est couramment utilisée. Elles
permettent de visualiser la conformation de la dose au volume cible et de localiser des régions de
surdosage ou sous-dosage.

D’un point de vue quantitatif, l’évaluation de l’histogramme dose volume (HDV) est envisagée. La
Figure 18 montre un HDV utilisé habituellement en clinique. Un HDV comprend généralement
l'ensemble des structures (volumes cibles et organes à risque) et des objectifs d'intérêt dans le plan de
radiothérapie, chaque courbe représente une structure spécifique. Chaque point d'une courbe
représente le pourcentage en volume de la structure associée recevant une dose au moins égale à la
dose lue en abscisse.

Figure 18 Un histogramme dose-volume (HDV) utilisé en clinique (Drzymala, et al., 1991)

Un HDV de la tumeur serait idéalement une fonction en échelon, avec 100% de la cible recevant
exactement la dose prescrite et 0% recevant plus. Un HDV idéal de la structure critique serait une
fonction avec 100% du volume de la structure recevant 0% de la dose prescrite. Un HDV capture les
informations volumétriques qui sont difficiles à percevoir à partir des courbes d'isodoses.

De nombreux indices sont proposés dans la littérature pour compléter ces évaluations. Nous allons
voir quatre types d’indices au sujet des différents volumes d’intérêt : indice d’homogénéité et l’indice
de conformation pour le volume cible (PTV), indice pour la protection de l’OAR, indice pour les tissus
sains (TS).

(ICRU 83, 2010) a recommandé un indice appelé « Homogeneity Index » (HI) pour déterminer
l’homogénéité globale de la dose absorbée par le PTV.

Il est défini comme indiqué ci-dessous:

Où :
représente la dose reçue par au moins du PTV.

Un très proche de 0 implique que la dose déposée dans le PTV est quasi homogène. L’objectif
d’optimisation est de faire décroitre l’indice vers la valeur idéale de 0.

31
Les indices suivants sont fondés sur les valeurs des isodoses de référence (IR), le volume de l’isodose
de référence et le volume d’intérêt (PTV, OAR). Ces différents volumes en question sont simplifiés en
2D et schématisé en Figure 19. Notons que (cercle vert), (cercle rouge) représentent le
volume du PTV, de l’OAR respectivement, (cercle bleu en pointillé) représente le volume de
l’isodose de référence, , représente le volume d’intérêt (PTV, OAR) couvert par l’isodose
de référence (IR).

Figure 19 Les volumes dans la définition des indices dosimétriques

Le Radiation Therapy Oncology Group (RTOG) a élaboré un index appelé « Conformation Index »
( ) pour déterminer la couverture du PTV par l’isodose de référence ( (ICRU 50, 1993) a
recommandé comme isodose de référence)) :

Il est défini comme indiqué ci-dessous :

La valeur idéale de est de 1. L’objectif d’optimisation est de faire croitre l’indice vers la valeur
idéale de 1.

De manière similaire, l’indice de protection de l’OAR est défini:

La valeur idéale de est de 1.

L’indice de protection des tissus sains est défini :

La valeur idéale de est de 0.

Dans la troisième partie du mémoire, nous allons évaluer notre plan de traitement par ces outils
présentés.

32
Deuxième partie

Calcul de dose

33
Chapitre I Dosimétrie ..................................................................................................................... 35
1 Transport de l’énergie par les photons et ionisation ...................................................................... 35
2 Le dépôt de l’énergie par les électrons .......................................................................................... 37
3 Terma, Kerma et dose absorbée .................................................................................................... 38
Chapitre II Les méthodes de calcul de dose : un état de l’art .......................................................... 41
1 Méthode Monte Carlo ................................................................................................................... 42
2 Méthodes basées sur la séparation des rayonnements primaires et diffusés ................................. 42
3 Méthodes de convolution / superposition ...................................................................................... 42
3.1 Méthode point Kernel ............................................................................................................ 43
3.2 Méthode Pencil Beam ........................................................................................................... 44
3.3 Méthode Collapsed cone convolution ................................................................................... 45
4 Méthodes de calcul par réseaux de neurones................................................................................. 46
Chapitre III Représentation analytique de la dose - Doséclair .......................................................... 48
1 Le maillage du fantôme ................................................................................................................. 49
2 Les modèles pré-calculés............................................................................................................... 50
3 Découpage du faisceau primaire en SBIMs .................................................................................. 50
4 Fonction de projection ................................................................................................................... 52
5 Calcul de la dose par somme pondérée de projections .................................................................. 53
6 L’algorithme de calcul de dose...................................................................................................... 58
Chapitre IV Contribution Doséclair .................................................................................................. 60
1 Les projections .............................................................................................................................. 60
1.1. Les repères de coordonnées ................................................................................................... 60
1.2. Le point origine du faisceau RCMI exprimé en coordonnée sphérique ................................ 61
1.3. Formulation d’une projection ................................................................................................ 62
2 Propagation axiale ......................................................................................................................... 64
2.1 Tenseur de la maille support.................................................................................................. 64
2.2 Simulation et résultat ............................................................................................................. 66
3 Propagation latérale ....................................................................................................................... 68
3.1 Tenseur de la maille voisine .................................................................................................. 68
3.2 Simulation et résultat ............................................................................................................. 70
4 Discussion ..................................................................................................................................... 72

34
 Chapitre I Dosimétrie
La traversée de la matière par un faisceau de particules aboutit à un dépôt d’énergie dans cette matière.
En radiothérapie, qui implique un faisceau de photons ou les électrons, le but de la dosimétrie est
d’évaluer quantitativement cette énergie absorbée afin de déterminer les effets des traitements sur les
tissus sains et les tissus tumoraux. Dans le cadre de cette thèse, nous nous sommes intéressés par des
faisceaux de photons de haute d’énergie.

1 Transport de l’énergie par les photons et ionisation

Les faisceaux de photons en radiothérapie sont des rayonnements indirectement ionisants qui met en
mouvement des particules chargées (électrons et positrons) dans le milieu; ces particules chargées, à
leur tour, vont interagir avec le milieu ce qui aboutit aux effets biologiques.

Pour les photons d'énergie inférieure à environ 10 eV, il n'y a généralement pas d'ionisation: le photon
ne possède pas suffisamment d'énergie pour arracher un électron aux atomes ou aux molécules de la
matière. Au-delà de 10 eV, l'énergie du photon est suffisante pour que l'ionisation soit possible.

Figure 20 (Ouabri, 2012) Les trois processus ionisants principaux des photons (a) photoélectrique (b) Compton (c)
création de paire

Les trois processus ionisants principaux illustrés en Figure 20 sont :

- l’effet photoélectrique (Figure 20(a)) :

L'énergie du photon incident est totalement communiquée à un électron périphérique de l'atome, qui
est éjecté de sa couche électronique.

- l’effet Compton (Figure 20(b)) :

L'énergie du photon incident est partiellement transférée à un électron périphérique de l'atome, qui est
éjecté de sa couche électronique. L'énergie restante est emportée par un photon diffusé.

- l’effet matérialisation (création de paire) (Figure 20(c)):

L'énergie du photon incident est totalement matérialisée sous la forme d’une paire électron-positron
sous l'effet du champ électrique du noyau atomique.

35
 Figure 21 Les modes d’interaction des photons en fonction de leur énergie (Johns, et al., 1983)

En radiothérapie, les interactions prépondérantes des photons de hautes énergies utilisées sont : effet
Compton entre 0.3-20 MeV et création de paires au-delà de 1.02MeV (Figure 21). Toutes ces
interactions conduisent à la production d’électrons de hautes énergies.

Lorsqu'un photon pénètre dans la matière, il n'est pas possible de prédire s'il va interagir et, s'il y a
interaction, où aura lieu celle-ci. On utilise la notion de probabilité pour tenir compte de cette
propriété: est la probabilité d'interaction d'un photon par unité de longueur de matériau traversé
(également appelé coefficient d'atténuation linéique). Le coefficient s'exprime en . Il ne dépend
que de l'énergie du photon et de la nature du matériau. Un photon a toujours la même probabilité
d'interagir dans la matière quelle que soit l'épaisseur de matériau traversé. Il a par conséquent une
probabilité non nulle de traverser la matière sans interagir.

Du coefficient d'atténuation on déduit le libre parcours moyen des photons, qui représente la
distance moyenne parcourue sans interaction par des photons d'une énergie donnée dans un matériau
donné. Ce libre parcours moyen , qui s’exprime en mètre, est donné par :

A titre d'exemple, le libre parcours moyen, dans l'eau, des photons de 3 MeV, est de 26 cm. Ceci
confirme le fort pouvoir de pénétration d'un faisceau de photons de haute énergie.

La probabilité d'interaction intervient dans la loi qui régit l'atténuation d'un faisceau de photons
incidents, sur un matériau d'épaisseur :

Le nombre de photons qui ont traversé sans interagir un matériau d'épaisseur (nombre de
photons transmis que l’on appelle la fluence primaire) est:

Les photons sont atténués exponentiellement avec la distance. Ainsi, le coefficient d'atténuation
linéique suffit à caractériser la fluence primaire en une position donnée.1

1
En toute rigueur, l’aspect probabiliste des interactions ne devrait pas permettre d’établir des valeurs scalaires
mais plutôt des densités. Toutefois compte tenu du nombre extrêmement importants de photons mis en jeu en
radiothérapie (plusieurs milliards au minimum), la variance statistiques est négligeable, indétectable avec

36
Bien qu'elle permette aisément le calcul de la densité d'interactions des photons primaires, la fluence
primaire ne peut suffire à caractériser la dose. L'énergie transférée par les photons n'est en effet pas
déposée localement, les électrons la transportant sur des longueurs non négligeables.

2 Le dépôt de l’énergie par les électrons

Dans le mécanisme de dépôt de l’énergie, les électrons secondaires mis en mouvement interagissent
avec les électrons périphériques des atomes du milieu, et par conséquent l’électron peut ioniser
plusieurs atomes sur son trajet.

Les principales interactions des électrons avec les atomes sont :

- la diffusion élastique : simple déviation de l'électron incident.

- la diffusion inélastique : qui provoque l'éjection d'un électron de l'atome, au prix d'une perte
d'énergie de l'électron incident.

- le rayonnement de freinage : l'électron est freiné et dévié en passant dans le voisinage d'un
noyau atomique. L'énergie perdue par l'électron est emportée par un photon, dit secondaire,
qui pourra interagir avec la matière plus loin. Ce processus bien que faible n’est pas
négligeable dans les calculs et explique que l’on observe un dépôt de dose non nul à
l’extérieur du faisceau, à une distance bien supérieure à la portée des électrons (voir ci-
dessous).

La trajectoire d’un électron est d’allure très chaotique et son parcours dans la matière est de longueur
bornée, comme on le voit sur la Figure 22. La distance moyenne de parcours d’un électron est
proportionnelle à son énergie initiale et inversement proportionnelle à la densité électronique du
matériau :

- plus l'énergie de l'électron est grande, c'est-à-dire plus il se déplace vite, plus son parcours sera
grand ;

- plus la densité électronique du milieu est grande, plus l'électron interagit, et plus son parcours
sera court.

Figure 22 Profondeur de pénétration (portée) d’un électron dans les tissus (Vuillez, 2009)

l’instrumentation classique en clinique, et donc on modélise par ces grandeurs scalaires qui en fait sont les
espérances.

37
Les électrons ont une portée bien plus courte que les photons. A titre d'illustration, la portée d'un
électron de 1 MeV est d'environ 0.4 cm dans l'eau, tandis que le libre parcours moyen d'un photon de
la même énergie est de 14 cm.

Autrement dit, les électrons sont vite absorbés et déposent toute leur énergie dans le tissu sous forme
d’ionisations; ceci permet d’atteindre des doses importantes nécessaires pour la stérilisation des
tumeurs.

3 Terma, Kerma et dose absorbée

Considérons au sein du matériau irradié une sphère élémentaire centrée sur de masse , illustrée
Figure 23. Pendant la durée de l’irradiation, un certain nombre de photons pénètrent dans la sphère, ils
transportent une énergie . Pendant la même durée, sortent de la même sphère des photons qui l’ont
traversé sans interaction. Soit l’énergie globale de ces photons qui sortent de la sphère. On pose :

représente l’énergie totale libérée. On appelle terma (Total Energy Released per Mass), noté , le
quotient par la masse de matière :

Parmi cette énergie totale libérée, une partie est transférée sous forme cinétique à des particules
chargées, noté . On appelle kerma (Kinetic energy released per mass unit), noté , le quotient
par la masse de matière :

L’unité du kerma est le gray, abrégé Gy, 1 Gy valant 1 J/kg.

Figure 23 Le kerma au sein de la sphère (Vuillez, 2009)

Les électrons mis en mouvement dans la sphère perdent leur énergie soit par collision soit en émettant
un rayonnement de freinage. On identifie la part du kerma due à ces deux phénomènes respectivement
comme le kerma de collision et le kerma radiatif. Notons que l’énergie transférée dans le kerma
collision n’est pas nécessairement absorbée dans la sphère car des électrons peuvent quitter celle-ci.

L’énergie cinétique cédée aux électrons à l’intérieur de la sphère va être absorbée en partie à
l’intérieur, en partie à l’extérieur de la sphère. Soit l’énergie absorbée à l’intérieur de la sphère,
aussi bien à partir des électrons mis en mouvement à l’intérieur ( et ) qu’à l’extérieur ( ),

38
illustré par la Figure 24. On appelle la dose absorbée, noté le quotient par la masse de
matière.

Dans l’eau, une dose absorbée de 1 Gy correspond à une élévation de température d’environ 2·10 −4 °C
(Vuillez, 2009).

Figure 24 La dose absorbée dans la sphère (Vuillez, 2009)

Dans la majorité des cas, seule une partie de l’énergie emportée par les électrons lors de l’interaction
en un point d’un rayonnement ionisant avec un milieu matériel est effectivement absorbée dans
l’élément de volume entourant le point ( et ), le reste de l’énergie transférée est absorbé
en dehors de l’élément de volume. A l’inverse, il absorbe tout ou partie de l’énergie d’électrons créés à
l’extérieur mais rentrant dans l’élément de volume ( ). On en déduit que par défaut le kerma de
collision et la dose absorbée sont différents. Mais sous certaines conditions, ces deux sources de
différences se compensent parfaitement et il y a égalité entre dose absorbée et kerma collision. Pour
établir cette relation, il faut que l’énergie emportée par les électrons secondaires sortant de l’élément
de volume considéré soit compensée par l’énergie déposée dans celui-ci par d’autres électrons
secondaires produit dans des éléments voisins. Cette situation caractérise l’équilibre électronique. Afin
que l’équilibre électronique soit réalisé, il faut que :

– Le volume de masse soit plongé dans un volume de matière homogène beaucoup plus
grand, dont les dimensions sont en rapport avec la portée des électrons.

– Le champ de rayonnement soit uniforme (ou puisse être considéré comme tel) sur .

Dans toute autre situation, l’équilibre électronique n’est pas assuré et il y a un écart entre le kerma
collision et la dose absorbée. Cela s’observe principalement aux interfaces entre différents matériaux
(non homogénéité sur ) et sur les bords du faisceau (non uniformité du champ de rayonnement sur
).

On définit comme le quotient de la dose absorbée rapportée au kerma collision (noté ), tel que :

Lorsqu’un photon de haute énergie pénètre le milieu, est maximum à la surface du matériau irradié,
car la fluence des photons est grande à la surface, puis décroit en profondeur. Ce processus va
engendrer une augmentation de la fluence des électrons secondaires, donc une augmentation de la dose
absorbée jusqu’à ce que la profondeur ou la dose atteint son maximum. À ce niveau, une
égalité de fluences d’entrée et de sortie est réalisée; fluence des photons (approximatif car la

39
décroissance en fluence n’est pas perceptible sur un petit volume) et fluence des électrons secondaires
; puis elles décroitront en même temps. Trois régions de sont définies :
: Augmentation de et diminution de .
: (équilibre électronique).
: Diminution de et parallèlement (équilibre électronique).

La variation du kerma collision et de la dose absorbée en fonction de la profondeur est représentée

par la courbe illustrée Figure 25.

Figure 25 Kerma collision et la dose absorbée en fonction de la profondeur dans un milieu (Gambini, et al., 2007)

Le kerma collision est proportionnel à la fluence qui a pour unité m-2 et au coefficient massique
d’absorption en énergie (unité m-2·kg-1), et à l’énergie du photon incident. Dans un milieu :

Par conséquent, pour un même faisceau de photon considéré projeté dans deux milieux différents, le
ratio de kerma collision entre le milieu et le milieu est :

Si l’équilibre électronique est vérifié (l’égalité entre le kerma collision et la dose) :

Compte tenu que les phénomènes de déséquilibre électronique sont similaires et vont dans le même
sens pour tous les matériaux, on peut étendre la formule au premier ordre et donc calculer la dose dans
un milieu donnée en connaissant pour le même faisceau de photon la dose dans un milieu de référence.

40
 Chapitre II Les méthodes de calcul de dose :
un état de l’art

En radiothérapie, il est indispensable d'avoir une connaissance précise de la dose délivrée au volume
cible et dans les organes critiques avoisinants. De nos jours, la dose est calculée en trois dimensions à
l’aide d’algorithmes implémentés sur des systèmes de planification de traitement (TPS) qui prennent
en compte les caractéristiques anatomiques des patients et les caractéristiques physiques et
géométriques des faisceaux.

En général, les algorithmes de calcul de dose modélisent les processus physiques de transfert d’énergie
et calculent la dose déposée par les faisceaux du traitement dans les voxels qui constituent le fantôme
numérique du patient. Selon (Ahnesjo, et al., 1999), la dose reçue par le patient est divisée en deux
catégories : premièrement, la dose venant des photons primaires qui n’ont pas eu d’interaction dans la
tête de l’accélérateur. Ces photons n’interagissent qu’avec les tissus du patient, qui provoque un dépôt
de dose par les électrons mis en mouvement et un dépôt (indirect) de dose par des photons diffusés
dans le fantôme. Deuxièmement, la dose venant des photons primaires qui ont interagit avec la tête de
l’accélérateur qui produit des électrons de contamination et les photons diffusés de la tête de
l’accélérateur. Les électrons de contamination déposent leur énergie à la peau du patient. Les photons
diffusés de la tête de l’accélérateur, à leur tour, déposent leur énergie dans le patient.

Ces algorithmes de calcul de dose pour des situations complexes (prise en compte des hétérogénéités,
des irrégularités de surface…) sont parfois très imprécis. Ils donnent une représentation approximative
de la répartition de dose dans le patient. Le niveau d’approximation sera différent selon le type
d’algorithme utilisé.

Le fait est qu’il est inévitable de faire des compromis entre la précision et la vitesse de calcul. À ce
jour, les algorithmes utilisés sont de plus en plus élaborés et visent à une amélioration de la précision
dans le calcul de la dose.

La précision des algorithmes de calcul de dose par le TPS est classiquement évaluée par la
comparaison entre les distributions de dose mesurées et calculées. Des grandeurs permettant d’évaluer
les écarts entre calculs et mesures sont définies. (Venselaar, et al., 2001) ont proposé de différencier
plusieurs régions de calcul en fonction du cas considéré pour distribuer des tolérances distinctes. Une
grandeur statistique limite de confiance basées sur les écarts de dose et l’écart type est définie pour
tenir compte d’une distribution de dose globale. Une autre méthode plus répandue, l’index gamma, a
été proposé par (Low, et al., 1998) pour effectuer une analyse globale d’une distribution de dose. Cet
index gamma combine différence de dose (en %) pour des zones à faible gradient et différence de
distance (en mm) pour des zones à fort gradient. L’utilisation courante en clinique prend des
différences de l'ordre de 2% et de 2 mm, ce qui est un bon compromis entre précision nécessaire
conservation d’un temps de calcul compatible avec une approche interactive.

L’index delta, introduit dans (Blanpain, et al., 2009), repose sur le concept de l’enveloppe delta. Cette
enveloppe transforme les tolérances en dose et en position de l’index gamma, en tolérances en dose
uniquement. La méthode est exposée en annexe F. Dans le cadre de cette thèse, l’index delta est utilisé
pour évaluer la précision du calcul de dose.

Les modèles de calcul de doses couramment utilisés sont ici succinctement décrits et commentés.

41
1 Méthode Monte Carlo
Les algorithmes Monte Carlo (MC) sont des méthodes stochastiques permettant de résoudre des
problèmes numériques pour lesquels on ne peut obtenir de formulation analytique. Dans le cas décrit
ici, les modèles physiques du processus de transport et de diffusion des électrons et photons sont en
fait des modèles statistiques d’interaction rayonnement-matière à l’échelle des particules (photons,
électron, atome). La formulation analytique complète d’un tel modèle est impossible, on lui applique
donc la méthode Monte Carlo (Andreo, 1991) (Salvat, et al., 1998). Mais par extension, en calcul de
dose, on appelle souvent « méthode Monte-Carlo » le duo modèle statistique et résolution
stochastique.

La méthode Monte Carlo a très souvent et très complètement été détaillée et commentée dans la
littérature. Elle est fondée sur la simulation du transport des particules depuis leurs production
jusqu’au dépôt de l’énergie dans la matière. Elle est la méthode actuelle la plus réaliste. Les
distributions de dose obtenues par la méthode Monte Carlo sont donc utilisées comme références.
C’est ce qui nous a amenés, au cours de ce travail de thèse, à utiliser des données pré-calculées par le
code Monte Carlo Penelope (Salvat, et al., 2006).

En revanche, pour obtenir des résultats précis avec une incertitude acceptable, il est indispensable de
simuler le transport d’un grand nombre de particules (au moins de l’ordre de 107). Ceci implique des
temps de calcul relativement longs qui peuvent atteindre plusieurs jours (Blanpain, 2009) (Menguy,
1996).

De nombreuses méthodes ont été développées pour accélérer le calcul de dose et permettre à terme
l’utilisation de cette approche dans la routine clinique. Parmi les principes les plus fréquents pour
permettre ces gains de temps de calcul, on peut noter les travaux sur la réduction de variance (Brualla,
et al., 2009), sur des approximations concernant les processus physiques (Habib, et al., 2010) ou
encore la parallélisation des calculs sur des clusters voire des cartes GPU (Love, et al., 2000). À ce
jour, le TPS Corvus v.5(NOMOS Corporation, Swickley, PA) a adopté la méthode Monte Carlo pour
le calcul de dose. Malgré tout, cette méthode reste encore peu adaptée à la dosimétrie numérique en
routine. Elles sont principalement réservées aux études théoriques, notamment pour tester la validité
des autres modèles dans des situations complexes, difficiles à réaliser sur le plan expérimental et
apparaissent dans les processus de validation clinique pour les cas complexes et critiques.

2 Méthodes basées sur la séparation des rayonnements primaires et

diffusés
Cette méthode a été longtemps utilisée dans les TPS au cours des années 90. Cette méthode a été
d'abord développée par (Clarkson, 1941), puis par (Cunningham, 1972) (Bjarngard, et al., 1989). Cette
méthode consiste à calculer séparément la « dose primaire » et la « dose diffusée ». La dose en un
point est la somme des contributions des composantes primaire et diffusée du champ de rayonnements.
La présentation détaillée se trouve dans (Menguy, 1996). Cette méthode est particulièrement adaptée à
la radiothérapie conformationnelle. De nos jours, elle est encore utilisée dans le TPS dédié du
Cyberknife, mais moins utilisée dans la RCMI.

3 Méthodes de convolution / superposition

42
 Figure 26 (Mackie, 2010)(a) Fluence du photon primaire (b) Noyau de convolution (c) Distribution de dose

Les méthodes de convolution / superposition, proposées au début des années 1980 (Boyer, et al., 1985)
(Mohan, et al., 1986) (Ahnesjo, 1989) (Mackie, et al., 1985), séparent les processus de transport et de
dépôt de l’énergie en deux phases : le transport de l’énergie par les photons primaires, et son dépôt par
les particules secondaires (électrons et photons). En réalité, les électrons issus des interactions
photoélectrique, Compton et création de paire déposent leur énergie dans les quelques millimètres
voire centimètres autour, mais les interactions Compton libèrent aussi des photons secondaires qui
peuvent parcourir de longue distance avant d’interagir de nouveau et de déposer toute ou partie de leur
énergie (Childress, et al., 2012). Le dépôt d’énergie par ces particules secondaires est également
intégré dans l’expression des Kernels. La dose absorbée totale en un point donné résulte de la
somme des distributions de dose calculées à partir de tous les points où l'énergie est libérée
(superposition), cette somme est réalisée par une convolution de ces deux phases, illustrée sur la
Figure 26 .

L’équation générale du calcul de dose s’écrit :

Où est le terma au point , donne la part de l’énergie déposée au point .

3.1 Méthode point Kernel

Figure 27: Noyaux de convolution : (a) noyau représentant la répartition de l’énergie libérée par des interactions de
photons primaires ayant lieu en un point unique ; (b) un noyau de type pencil beam, représentant la dose déposée par
un faisceau de section infinitésimale.

La méthode du point kernel consiste en deux phases successives. La première phase correspond à
calculer l’énergie totale transférée par tous les photons primaires dans une unité de masse (terma, vue

43
en section 3 Chapitre I de la première partie du mémoire). La deuxième phase correspond à calculer un
modèle de dépôt de cette énergie autour d’un site d’interaction primaire.

Puisque le déplacement des électrons et le dépôt d'énergie est aléatoire, le modèle de dépôt de cette
énergie est évalué en utilisant des simulations Monte Carlo. Ce modèle est appelé noyau de
convolution (kernel : point spread functions (Ahnesjo, et al., 1987)) qui décrit la distribution de
l’énergie déposée dans un milieu infini homogène autour d’un site d’interaction primaire. Il peut être
observé sur la Figure 27(a).

Il faut noter que les noyaux de convolution sont généralement calculés dans un fantôme d'eau, alors
que les patients sont composés de milieux divers comme des tissus musculaires, de l’os, du poumon,
etc. Dans la pratique, on considère que les photons de haute énergie interagissent d’une manière
similaire dans le poumon, les tissus musculaires, les tissus adipeux et même l’os que dans l'eau.
Cependant, les photons de basse énergie, grâce à l’effet photoélectrique, sont plus facilement absorbés
dans un milieu dense comme l’os par rapport à l’eau. Des matériaux avec un grand numéro atomique
(comme le titane dans une prothèse) peuvent avoir aussi des proportions d’interactions très différentes
par rapport à celles calculées dans l'eau.

Pour corriger ce biais, l’adaptation aux cas hétérogènes a souvent été réalisée sur les deux phases du
calcul:

Tout d’abord, le terma qui est proportionnel à la fluence primaire doit être calculé en prenant en
compte les coefficients d'atténuation des différents matériaux traversés. Ensuite, le noyau de
convolution doit être dilaté ou compressé en fonction de la variation de la densité électronique. Ces
compressions et dilatations sont communément dénommées scaling (O’Connor, 1984). L’équation
générale du calcul de dose s’écrit :

Où est la densité électronique du milieu où le noyau a été généré, est la densité

électronique du milieu où le noyau doit être adapté.

La méthode de convolution est coûteuse en temps de calcul, puisqu'elle nécessite deux boucles
imbriquées sur l'ensemble des voxels du fantôme. A titre d’exemple, pour 106 voxels (fantôme de
30×30×30 cm3 avec une résolution de 3 mm), il faut 1012 calculs, ces calculs durent donc plusieurs
minutes pour un fantôme 3D. Pour accélérer le calcul, différents méthodes ont été envisagées. La
première est le calcul de la convolution par transformation de Fourier (Mohan, et al., 1987). Mais cette
méthode ne peut pas être adaptée à l’hétérogénéité du fantôme car elle utilise un noyau constant. La
seconde approche accélératrice consiste à pré intégrer des noyaux dans certaines directions. Ces
méthodes sont les plus efficaces et sont souvent implémentées dans les TPS commercialisés en
clinique, la plus répandue étant la méthode Pencil Beam décrite à la section suivante.

3.2 Méthode Pencil Beam

Une pré-intégration de tous les noyaux dans le sens de la profondeur du faisceau de section
infinitésimale permet d’obtenir la dose absorbée due à ce dernier. Un tel noyau est montré sur la
Figure 27 (b). Cette approche dite pencil beam est utilisée pour les faisceaux d’électrons (Kooy, et al.,
1989) et les faisceaux de photons (Ahnesjo, et al., 1992).

Pour calculer la dose déposée par un faisceau plus large, tous ces noyaux sont alors mis côte à côte.
On réalise ainsi une convolution entre le noyau pencil beam et le terma du faisceau dans le fantôme.
Cette méthode donne d’excellents résultats dans un milieu homogène, l’adaptation aux cas

44
hétérogènes a été réalisée en faisant un scaling sur le noyau dans un milieu de référence (souvent de
l’eau) selon la densité électronique du milieu courant. Cette adaptation est limitée par un scaling selon
la densité électronique le long de l’axe du faisceau, ce qui fait que les diffusions latérales sont
généralement négligées. En outre, cette méthode ne permet pas de calculer les variations de dose aux
interfaces issues du déséquilibre électronique local. Par conséquent, les résultats obtenus peuvent être
très imprécis, des erreurs importantes ont été relevées (Knoos, et al., 1995).

Cependant, des améliorations ont récemment été apportées à la méthode du pencil beam (Ulmer, et al.,
2005) (Tillikainen, et al., 2008). Le problème de diffusion latérale est résolu par un scaling latéral des
noyaux, et des modèles complémentaires ont été ajoutés aux interfaces pour compenser les
déséquilibres électroniques induits par celles-ci.

Grâce à ces deux améliorations, l'approche du pencil beam devient intéressante pour les applications
cliniques avec un temps de calcul assez réduits. A titre exemple, le temps de calcul peut aller de dix
secondes à une minute pour des tailles de champ allant de 4×4 cm² à 6×6 cm².

La précision a par ailleurs été fortement augmentée par les dernières améliorations (Tillikainen, et al.,
2008) mais des erreurs importantes (jusqu'à 8% dans le cas de faisceaux étroits) peuvent toutefois être
relevées sur des géométries complexes.

3.3 Méthode Collapsed cone convolution

La méthode dite ''Collapsed cone convolution'', introduite par (Ahnesjo, 1989) offre un des meilleurs
compromis temps/précision.

Figure 28 Un voxel terma (cube bleu) émet un cône d'énergie. La dose est déposée seulement dans les voxels qui sont
sur l’axe de ce cône (flèche)

Pour un voxel avec un terma donné, la méthode considère que le transport d’énergie se fait selon des
cônes dans les différentes directions et partant du point central (le point où le terma a été
préalablement calculé) (voir Figure 28). Ensuite, elle suppose que toute l’énergie qui est propagée
dans un cône est transportée, atténuée et déposée selon une loi exponentielle sur l’axe de ce cône.
Cette hypothèse permet de simultanément distribuer la dose le long d’un ensemble de raies discrètes
qui émerge de chaque point terma et d’accumuler le terma au fur et à mesure que l’on avance dans la
direction considérée (voir Figure 29). Cette simultanéité est à la base des performances en temps de
calcul de l’algorithme.

45
 Figure 29 La propagation de l’énergie d'un voxel terma (cube) le long des cônes

La méthode qui en résulte est aujourd'hui l'une des plus efficaces, tant elle offre une précision
acceptable (rarement plus de 5% d’erreur (Vanderstraeten, et al., 2006) (Fogliata, et al., 2007)), tout en
prenant un temps limité (de l’ordre de la minute). Mais ce temps est toutefois encore trop élevé pour
envisager l’optimisation itérative des plans de traitement.

4 Méthodes de calcul par réseaux de neurones

Ces dernières années, différents auteurs ont utilisé des réseaux de neurones, tirant profit de leurs
propriétés d’approximateurs universels (Dreyfus, et al., 2002) pour calculer rapidement et précisément
la dose. Il est montré qu’un réseau de neurones est capable « d’apprendre » la dose en milieu
homogène sur l’axe du faisceau (Wu, et al., 2000) et la dose sur un plan 2D dans des volumes
homogènes (Blake, 2004). L’approche la plus élaborée à ce jour, nommée NEURAD (Mathieu, et al.,
2005) (Bahi, et al., 2006) est capable de calculer la dose absorbée en milieu hétérogène par un faisceau
large à partir de doses en milieux homogènes préalablement apprises par réseaux de neurones. Cette
méthode est fondée sur les deux hypothèses suivantes :

1. Dans une sous-partie homogène d'un fantôme hétérogène, il est supposé que la dose peut être
obtenue à partir d'un modèle pré-calculé en milieu homogène, et appris par un réseau de
neurones.

2. L'équilibre électronique est supposé réalisé au voisinage des interfaces (des interfaces
traversées par tout ou partie du faisceau primaire). La conséquence de cette hypothèse est que
la dose ne subit pas de fortes variations près de ces interfaces.

Cette méthode calcule la dose dans le cas d’une interface orthogonale à l’axe du faisceau. Elle
considère qu’il y a égalité entre la dose au dernier point du premier matériau et la dose au premier
point du second. Ils utilisent les courbes de dose, obtenues via leur réseau de neurones, dans le premier
et le second matériau. La dose avant l’interface est considérée égale à la dose en milieu homogène. La
dose après l’interface l’est également, mais elle est prise à partir de la profondeur où elle est égale à la
dose au dernier point avant l’interface. De cette façon, la courbe de dose à l’interface est bien
continue.

Les résultats de cette méthode sont très corrects dans le cas d'un faisceau large. En revanche, le cas des
faisceaux étroits n'est pas correctement traité (en fait la méthode exploite le théorème de Fano2 qui

2
(AAPM, 2004) Dans un milieu suffisamment étendu, possiblement hétérogène en densité, mais homogène en
composition, si la fluence en photons primaires est la même partout, alors la dose est constante. Or, au centre
d’un faisceau large, la fluence peut être considérée constante dans une zone suffisamment étendu autour de l’axe
d’un faisceau. On la considère constante avec la profondeur avant et après l’interface, son atténuation étant très
lente. Les conditions du théorème sont réunies, ce qui explique la faible variation de dose à l’interface. En

46
n’est pas valide par hypothèse dans le cas de faisceaux étroits). Une méthode de traitement des
interfaces parallèles à l'axe du faisceau, et à l'intérieur de celui-ci, est également donnée (Mathieu, et
al., 2005) (Bahi, et al., 2006). Néanmoins, la question des interfaces obliques, ainsi que celle des
interfaces en dehors du faisceau, n'est pas traitée. Le temps de calcul sur des grilles 3D est supérieur à
la minute.

Malgré les limitations énoncées ci-dessus, la méthode a contribué à l'idée d'utiliser, dans une sous-
partie homogène du fantôme, une distribution de dose pré-calculée dans un fantôme homogène
constitué du même matériau. La question majeure est de gérer le déséquilibre électronique près des
interfaces dans le cas de petits faisceaux. Mais un autre point important est également à considérer
pour avoir les meilleures performances.

Les méthodes présentées jusqu’à présent calculent des distributions de dose complètes. En d'autres
termes, si l'on ne souhaite calculer la dose qu'à une position unique, il est nécessaire de faire de
nombreux calculs intermédiaires sur d'autres points, voire de faire un calcul de dose complet sur une
importante partie du fantôme.

Or, obtenir une distribution de dose complète sur un fantôme complet n'est pas toujours nécessaire. Par
exemple, si la dose doit être connue de façon complète et précise sur la tumeur et les organes à risques,
on peut en revanche diminuer la densité de points de contrôle sur les zones moins sensibles. Un autre
exemple est l'optimisation du plan de traitement, pour laquelle il peut être suffisant, surtout dans les
étapes préliminaires de placement et d'orientation des faisceaux, de ne disposer que de la dose en
certains points. Par exemple, quelques points de contrôle peuvent suffire à vérifier si l'orientation d'un
faisceau est correcte, ou s'il convient de la modifier.

La méthode Doséclair, décrite ci-après et qui a été totalement reformulée au cours de la première
partie de cette thèse a l'avantage de fournir un résultat en un temps à peu près proportionnel au nombre
de points où la dose est demandée. Ceci est un avantage majeur par rapport aux méthodes existantes.

revanche, dans le cas d’un faisceau étroit, la fluence primaire n’est pas constante sur une assez grande largeur.
Les flux électroniques latéraux ne sont pas équilibrés.

47
 Chapitre III Représentation analytique de la
dose - Doséclair
La méthode Doséclair a été développée dans le cadre d’une thèse (Blanpain, 2009) qui s’est déroulée
dans le Laboratoire Information Modèles et Apprentissage du CEA LIST. Cette méthode modélise le
dépôt de dose d’une façon déterministe et calcule une distribution de dose dans un fantôme 3D
segmenté en mailles homogènes.

La méthode réalise ce calcul en segmentant le fantôme du patient en mailles homogènes, et en

associant aux mailles des projections vers des distributions pré-calculées en milieux homogènes, ainsi
que des pondérations gérant les hétérogénéités.

Le calcul de dose est réalisé en deux étapes. La première étape concerne les mailles : les projections et
pondérations y sont paramétrées en fonction de critères physiques et géométriques. La seconde étape
concerne les voxels : la dose y est calculée en évaluant les fonctions préalablement associées à leur
maille.

Cette méthode est très rapide, notamment quand le nombre de points d'intérêt est limité (quelques
centaines), les résultats étant dans ce cas obtenus en moins d'une seconde. Avec de telles
performances, la planification automatique des traitements de radiothérapie devient parfaitement
envisageable.

Un schéma de synthèse de différentes étapes de la méthode est présenté ci-dessous (Figure 30).
L’algorithme commence par segmenter le fantôme hétérogène en mailles homogènes (section 1) suivi
par un pré-calcul de dose déposée dans chaque modèle homogène prédéfini (section 2). Ensuite, le
faisceau est découpé en SBIMs (section 3) pendant sa propagation dans le fantôme suivi par des
calculs de fonction de projection (section 4) et pondération (section 5). En fin, la dose est calculée
pour les points d’intérêt. Le temps de calcul est analysé en section 6.

Maillage du fantôme

Préparation modèle homogène

Découpage en SBIMs

Projection Pour chaque SBIM Pondération

Calcul de dose par somme pondérée

Figure 30 Les étapes de la méthode Doséclair

48
1 Le maillage du fantôme
Par rapport à des méthodes comme Monte-Carlo, les approches telles que pencil beam et NEURAD
utilisent un changement d’échelle sur le faisceau (on passe des doses des particules individuelles à
celles des faisceaux complets) pour accélérer les calculs. Pour améliorer ce dernier, Doséclair utilise
en plus un changement d’échelle sur le fantôme hétérogène de la même nature (des voxels aux
mailles).

En effet, un corps humain n'est pas constitué de voxels de différents matériaux distribués de façon
anarchique, mais plutôt de tissus, d'os et d'organes. Il est donc naturel de chercher à regrouper les
voxels du fantôme en des ensembles homogènes -- les mailles --, et de chercher à les utiliser
avantageusement pour le calcul de dose.

Le regroupement de voxels adjacents en groupes homogènes est appelé segmentation. De nombreux

algorithmes existent pour réaliser cette tâche (Horaud, et al., 1995). Ils se regroupent en différentes
catégories, selon qu'ils considèrent principalement les régions homogènes ou leurs contours.

La représentation numérique du fantôme maillé est un des points clés de la méthode proposée dans ce
mémoire. Pour limiter la complexité de la représentation du fantôme maillé, la structure d’une maille
doit être décrite par un faible nombre de paramètres. Les mailles peuvent par exemple être des sphères,
des cubes, ou des parallélépipèdes rectangles.

Le prototype logiciel développé (Doséclair) dans le cadre de ce travail s’appuie sur des maillages
réguliers en parallélépipèdes rectangles, illustré en Figure 31.

Figure 31 Un maillage régulier de parallélépipèdes rectangles (illustration en 2D)

La régularité du maillage permet de le coder sous la forme d'un tableau 3D dont chaque case
correspond à une maille. Il est ainsi très rapide d'identifier la maille à laquelle appartient un point du
fantôme. La conception informatique de ces traitements se trouve donc simplifiée. Toutefois pour
l’étude théorique, nous nous autoriserons un maillage non régulier plus générique.

Ce maillage a l'inconvénient de ne pas être adaptable de façon simple à des anatomies présentant des
formes arrondies, pourtant très présentes dans l'anatomie humaine. Pour ce faire, des formes arrondies
peuvent être représentées soit par de nombreuses mailles rectangulaires de petites dimensions, et donc
au détriment de la vitesse des calculs, soit par peu de mailles rectangulaires, et donc au détriment de la
précision géométrique. Nous verrons plus tard que le travail de reformulation de la méthode fait dans
le prochain chapitre ouvrira aussi de nouvelles perspectives.

Ces questions n'ont pas été traitées au cours de ce travail. Elles devront l'être au cours de l'étude
complémentaire qui adaptera la méthode aux conditions cliniques réelles.

49
Le maillage parallélépipédique choisi convient cependant parfaitement au démonstrateur logiciel,
celui-ci ayant pour but de montrer qu'un maillage peut être utilisé avantageusement pour le calcul de la
dose.

2 Les modèles pré-calculés

La méthode Doséclair exploite des distributions pré-calculées de dose dans les fantômes homogènes
qui sont constitués des différents matériaux que le faisceau est susceptible de rencontrer dans un
patient. À titre exemple, dans ce mémoire trois distributions de dose dans un fantôme homogène sont
utilisés : dans de l’eau, dans de l’os, dans du poumon.

Le faisceau utilisé dans ces pré-calculs en milieux homogènes est le même que l'on simule ensuite
dans un fantôme hétérogène quelconque (mêmes dimensions, et même spectre en énergie). À titre
d'exemple, le faisceau utilisé dans les tests proposés au chapitre IV est un faisceau élémentaire de
photons de RCMI (section de1×1 cm2), mono énergétique de 5 MeV. En revanche, la seule différence
est que la fluence du faisceau utilisé en milieux homogènes est unitaire, alors que ce n’est pas
forcément le cas pour la fluence du faisceau utilisé dans un fantôme hétérogène. La dose finale
calculée dans le fantôme hétérogène est multipliée par la vraie fluence du faisceau pour le fantôme
hétérogène.

Issue de ces pré-calculs, un ensemble de modèles sont obtenue. , au point , la dose

déposée par le faisceau avec une fluence unitaire dans un fantôme homogène composé du
matériau . Nous allons calculer la distribution de dose du fantôme hétérogène à l’aide de ces
modèles pré-calculés.

En pratique, on utilise des fantômes homogènes suffisamment larges pour que la dose ne soit pas
significative à leurs bords. Les pré-calculs sont réalisés à l’aide du code Monte Carlo Penelope
(Salvat, et al., 2006).

3 Découpage du faisceau primaire en SBIMs

Nous avons vu que les méthodes de convolution /superposition (point kernel, collapsed cone
convolution) calculent, en tout point, l'énergie primaire libérée par les photons primaires (terma). Cette
énergie est ensuite répartie autour des points selon des noyaux de convolution, pour donner la dose. La
dose totale délivrée par le faisceau est alors calculée comme la somme des contributions de chacun des
voxels où passe le faisceau.

Le calcul réalisé par ces méthodes utilise un fractionnement du faisceau basé sur les voxels. Or, une
méthode très rapide ne peut pas être obtenue en utilisant une granularité aussi fine. Nous devons donc
procéder autrement.

C'est ici qu'intervient le maillage du fantôme : grâce à lui, au lieu de découper le faisceau sur la base
des voxels, nous proposons de le faire à partir des mailles. Plus précisément, nous considérons, pour
une maille donnée, des sous-parties (une ou plusieurs) du faisceau qui passent dans les mailles. On
appelle SBIM (Sub Beam In a Mesh) une sous-partie du faisceau primaire passant dans une maille
donnée. Ce concept est présenté à travers l’exemple de la Figure 32, qui montre le découpage du
faisceau primaire en cinq SBIMs. On remarque en particulier sur cet exemple qu’il peut y avoir un ou
plusieurs SBIMs dans une même maille (par exemple, deux SBIMs dans la maille 4 en bleu à droite).

50
 Figure 32 Découpage d’un faisceau en SBIMs, lors de son passage dans un fantôme hétérogène à quatre mailles

Avec la même analogie que la méthode de convolution/superposition, la dose totale déposée par le
faisceau est la somme des contributions de chacun de ces SBIMs.

La méthode Doséclair impose que le profil de fluence de chaque SBIM doit être équilibré, c’est-à-dire
constant ou d’une très faible variation sur sa section. Le profil de fluence d’un SBIM peut se trouver
déséquilibré à la suite d’un chevauchement de matériaux différents ou la traversée oblique d’une
interface. Dans le premier cas, le profil de fluence est hétérogène composé de plusieurs parties
homogènes (des plateaux avec un SBIM par plateau). C’est le cas des SBIM 4 et 5 sur la Figure 32 car
le faisceau a traversé les mailles 2 et 3 composées de matériaux différents. Dans le deuxième cas, le
profil de fluence est également déséquilibré, mais cette fois-ci en continu. La solution consiste à
découper le faisceau en SBIMs suffisamment étroits pour que les variations de leur fluence ne soient
pas significatives sur la section. Cette approximation sur la fluence est illustrée sur la Figure 33.

Figure 33 Découpage d’un faisceau en deux SBIMs (à peu près homogènes en fluence) lors d’une traversée d’interface
oblique).

Le calcul de SBIMs via le maillage est fait en deux temps : d’abord le découpage lors de la
propagation du faisceau primaire dans le fantôme, puis une phase d’agrégation des SBIMs inutilement
séparés.

Le découpage d’un faisceau en SBIMs se fait dans les deux cas suivants :

1. le faisceau pénètre dans plusieurs mailles à la fois, illustré sur la Figure 34 où le faisceau est
séparé en deux à son entrée dans le fantôme.
2. le faisceau traverse une interface oblique et il est assez large pour qu’un découpage en largeur
soit nécessaire compte tenu de l’incidence (voir Figure 33).

51
Figure 34 Découpage d’un faisceau traversant obliquement un fantôme maillé (la numérotation des SBIMs ne
représente pas nécessairement leur ordre de création)

Ce processus est répété chaque fois qu’un SBIM se présente à l’entrée d’une maille : soit il est
transmis, soit il est découpé.

La mise en œuvre de cette propagation par le logiciel nécessite la gestion des contours des SBIMs
(leur section). Lorsque le faisceau et les mailles sont délimités par des arêtes rectilignes, ces sections
peuvent être modélisées par des polygones.

En se propageant dans le fantôme, le faisceau primaire peut être découpé un grand nombre de fois. Il
peut en résulter un grand nombre de SBIMs sur une même maille, avec en particulier des SBIMs dont
la séparation n'a pas réellement lieu d'être, illustré en Figure 35. Le SBIM 1 et SBIM 2 sont créés
pendant la propagation du faisceau, mais cette séparation n'est pourtant pas nécessaire à cette
profondeur, car ces deux SBIMs ont exactement la même fluence (la partie haute et la partie basse du
faisceau ont traversé autant d'os et de poumon l'une que l'autre).

Ces découpages inutiles ont en effet comme conséquence une augmentation du nombre de modèles à
évaluer lors du calcul de la dose sur la maille concernée. Cette augmentation allonge la durée du
calcul, sans améliorer sa précision. La solution est de fusionner les SBIMs adjacents dont la séparation
n'a pas lieu d'être à chaque fois qu’une maille est traitée.

Figure 35 Illustration d’une agrégation possible des deux SBIMs

4 Fonction de projection
Le fantôme hétérogène est segmenté en mailles homogènes. Grâce à ces mailles homogènes, le
faisceau primaire est segmenté en SBIMs. Les modèles pré-calculés nous donnent la dose déposée en
chaque point du fantôme homogène.

52
L’idée de base de la méthode Doséclair est de trouver un SBIM virtuel dans le fantôme en milieu
homogène associé, pour chaque SBIM créé dans le fantôme patient. La distribution de dose dans ce
SBIM virtuel est équivalente à celle du SBIM réel dans le fantôme hétérogène. En utilisant une
fonction de projection , on peut établir une correspondance entre le SBIM réel et le SBIM virtuel.

Un exemple est illustré en Figure 36, un SBIM dans le fantôme hétérogène à gauche est projeté dans
son modèle homogène à droite. La dose déposée en point (resp. ) est équivalent à la dose déposée
en (resp. ). Cette projection doit être réglée de sorte que la fluence primaire sur le
SBIM ait autant d’atténuation que la fluence primaire de son projeté dans le fantôme homogène.
Ramené à un point, cela signifie que la profondeur du point (resp. ) dans le modèle
homogène est telle que la fluence locale en ce point (resp. ) dans ce modèle
homogène est la même que la fluence locale en le point (resp. ) dans le fantôme hétérogène. Cette
propriété doit en particulier être vérifiée au centre du SBIM.

Dans la suite de ce mémoire, l’exposant est utilisé pour désigner les SBIMs résultant du découpage
d’un faisceau. La projection associée au SBIM , représenté par , est une fonction analytique.
La dose obtenue au moyen de cette projection est notée :

est la distribution de dose dans le fantôme homogène de matériau . Le calcul de

revient à projeter le point du fantôme patient en un point du fantôme
homogène, puis à évaluer le modèle en ce point.

Nous détaillerons le calcul de dans le Chapitre IV.

Figure 36 La fonction de projection d’un SBIM associe, à chaque position du fantôme hétérogène, une position du
fantôme homogène du même matériau que la maille support du SBIM

5 Calcul de la dose par somme pondérée de projections

La question est maintenant de savoir comment utiliser les projections associées aux SBIMs pour
calculer la dose totale. Si on suppose simplement que la dose, dans un SBIM du fantôme maillé, est
égale à la dose dans la projection de ce SBIM sur le modèle homogène pré-calculé, alors pour un point

53
 situé à l’intérieur du SBIM , la dose serait . Un tel mécanisme impliquerait que le
dépôt total de la dose est strictement local, ce qui donnerait une dose nulle en dehors du faisceau et
une discontinuité de dose aux frontières entre les SBIMs.

Figure 37 Un faisceau composé de 4 SBIMs dans un fantôme maillé

En réalité, la diffusion des électrons disperse l’énergie libérée lors des interactions primaires. Les
SBIMs ont également de l’influence sur la dose dans leur voisinage. En conséquence, la dose absorbée
en un point donné peut provenir d'interactions primaires ayant eu lieu dans différents SBIMs. Pour
illustrer cette affirmation par les méthodes classiques, il suffit d’envisager la dose d’un SBIM en la
calculant par la méthode à noyau. Par définition du SBIM, le terma sera non nul uniquement dans le
volume défini par le SBIM. Mais en convoluant par le noyau, on arrive immédiatement à la conclusion
qu’un SBIM influence la dose en dehors de sa localisation stricte, comme tous les SBIM adjacents.

Sur l’exemple de la Figure 37, un point situé sur l’interface séparant les SBIMs 1et 3 peut recevoir de
la dose venant de ces deux SBIMs. La dose en un point donné devra par conséquent être calculée à
partir des projections de plusieurs SBIMs, c’est-à-dire de plusieures contributions . Pour ce
faire, une fonction de pondération associée à chaque SBIM représente la part de
à la dose absorbée totale en . Dans le cas de notre exemple du point à l’interface entre
les SBIMs 1 et 3, les travaux antérieurs ont montré qu’une bonne pondération était et
. Nous allons voir maintenant comment systématiser leur calcul.

D’une façon générale, la dose totale, en un point quelconque, sera calculée de la manière suivante :

Où est la fluence du faisceau à son entrée dans le fantôme maillé.

La pondération dépend de la position et de la forme du SBIM. La dose totale est un

mélange des différentes . Par conséquent, la somme des pondérations , associées
aux différents SBIMs, doit être égale à 1 en tout point du fantôme :

54
Dans la pratique, comme nous le verrons, il peut y avoir plusieurs concepts physiques simultanés qui
influencent et permettent donc de calculer ces coefficients. Il n’est toutefois pas trivial de considérer
l’ensemble des phénomènes simultanément tout en garantissant la normalisation à 1. Nous allons donc
plutôt définir des « lois d’influence » selon des fonctions , et les coefficients en seront déduits
par la normalisation :

Le fantôme est composé des mailles homogènes. Pour un SBIM donné, on va appeler maille support la
maille où est localisé le SBIM et mailles voisines les autres mailles, suffisamment proches pour que le
SBIM considéré contribue significativement à la dose déposée (et donc au calcul).

Dans la maille support, le calcul de à partir des et des SBIMs. Pour un

SBIM , est calculé à partir de la fonction de pondération latérale et pondération axiale.

Sur l’exemple de la
Figure 38, à l’entrée du fantôme, le faisceau primaire est segmenté en deux SBIM, ayant chacun une
fluence uniforme (Figure 38(a)). Sur la
Figure 38 (b), les doses (en bleu) et (en rouge) sont deux demi-profils associés
respectivement aux SBIMs 1 et 2. Elles sont obtenues par leur fonction de projection :

La dose en noire est la référence calculée par Penelope. Nous pouvons constater que, sur cet exemple,
la projection seule donne une dose correcte sur les côtés mais une discontinuité au centre. Ceci est dû
au fait que les SBIMs sont considérés n’avoir qu’une action locale. En réalité, un SBIM dépose de la
dose en « dehors » de lui-même3. Sur la
Figure 38 (c), les projections de la dose sont étendues à toute la largeur. Sur la
Figure 38 (d), une pondération latérale de type gaussienne (avant normalisation) est mise en place.
Grace à cette gaussienne, chacun des deux SBIMs contribue majoritairement à la dose dans sa zone.
La discontinuité de la
Figure 38 (b) est remplacée par une transition douce de la dose entre ces deux zones, qui s’explique
par les diffusions électroniques latérales. Le résultat est montré en
Figure 38 (e), la somme des doses issues des projections des SBIMs pondérées par des gaussiennes
donne un bon accord avec la dose de référence calculée par Penelope.

3
Pour illustrer ceci avec des modèles connus, reprenons les méthodes par convolution. En fait,
géométriquement, c’est le terma qui coïncide exactement avec le SBIM. Dès que l’on convolue ce terma par un
noyau, la zone où la dose n’est pas nulle inclut mais est plus grande que la zone de terma non nul (donc le
volume du SBIM).

55
 Figure 38 Illustration de la pondération latérale

Sur l’exemple de la Figure 39, le faisceau primaire traverse deux mailles de compositions différentes.
Il passe d’abord dans une maille d’eau puis une maille d’os (Figure 39 (a)). Un SBIM est créé sur
chacune de ces deux mailles. Les projections associées aux SBIM 1 et SBIM 2 se font respectivement
vers les modèles eau et os :

On voit sur la Figure 39 (b), le rendement en profondeur dans l’eau en bleu et dans l’os en rouge. La
courbe noire est la dose calculée par les modèles de projection sans pondération. Rappelons que les
projections sont calculées de sorte que, pour un point de l’interface, la fluence primaire soit la même
au point du modèle d’eau, et au point du modèle d’os. est la profondeur à
laquelle la fluence primaire du modèle d’os est égale à la fluence primaire du modèle d’eau à la
profondeur (interface). La dose noire dans la partie droite est obtenue par simple décalage sur la
dose dans l'os.

Sur la Figure 39 (c), le résultat ainsi obtenu est comparé avec la dose de référence Penelope. Notre
dose calculée par une simple projection est discontinue à l’interface sachant qu’en réalité, après
l’interface, il y a un retour progressif à l’équilibre électronique selon la courbe de référence. Cette
zone de discontinuité soulignée par un cercle rouge est zoomée en Figure 39 (d). Sur la Figure 39 (e),
les deux courbes de dose sont étendues dans la maille mitoyenne. Deux pondérations axiales en
sigmoïdes (de la forme ) sont appliquées sur leurs courbes de dose associées pour modéliser la
transition au voisinage de l’interface. Les pondérations sont paramétrées pour passer de 0 (resp. 1)
quelques millimètres avant l’interface, à 1 (resp. 0) un centimètre après cette interface. Sur la Figure
39 (f), on voit que la transition réalisée (en noir) est assez proche de la référence (en vert). À noter

56
qu’ici la pondération a été directement faite sur les résultats de Monte-Carlo pré-calculés en milieu
homogène, sans aucun lissage ou apprentissage.

Figure 39 Illustration de la pondération axiale

En pratique, la pondération de forme gaussienne et la forme sigmoïde nécessitent une évaluation de la

fonction exponentielle, qui est coûteuse en temps de calcul. Pour une raison de performance, les
pondérations, dans (Blanpain, 2009) , sont modélisées par une fonction ayant un plateau de valeur 1 et

57
une croissance (resp. décroissance) de la forme à l’extrémité du plateau, illustré par Figure 40.
Notons que est la distance à l’extrémité du plateau le plus près, et est un paramètre permettant
d’ajuster la pente de la courbe.

Figure 40 La pondération implémentée en Doséclair

Pour un cas complet, la pondération d’un point du SBIM s’écrit :

Où est la pondération axiale, est la pondération latérale.

Rappelons que les coefficients s’écrit :

Une fois que les modèles complets ( et ) des SBIM sont établis, ils sont étendus
aux mailles voisines où il n’y a pas de SBIM directement créé. Pour une maille voisine composée du
même matériau que la maille support du SBIM en question, l’extension est directe puisque l’on garde
la même fonction de projection et les mêmes fonctions de pondération. Dans le cas contraire, un
scaling latéral est réalisé dans la projection pour prendre en compte la différence de
décroissance latérale de la dose due au changement de matériau. Cette technique de scaling sera
détaillée dans le prochain chapitre.

6 L’algorithme de calcul de dose

Avant d’utiliser la méthode pour calculer la dose due à un faisceau dans un fantôme hétérogène maillé,
il est nécessaire de précalculer les modèles de dépôt de dose de ce même faisceau en milieux
homogènes. Ces modèles sont les fonctions décrites dans la section 2. Une fois que ces
modèles sont disponibles, on peut passer au calcul de dose dans ce fantôme.

L’algorithme de calcul de dose est illustré sur un exemple, Figure 41, et suit les deux grandes étapes
suivantes :

- Préparation des modèles dans les mailles :

Le faisceau est propagé de maille à maille, les SBIMs sont créés (la Figure 41 (b)). Pour
chaque SBIM , une fonction de projection et de pondération sont définies,

58
 illustré par la Figure 41 (c). Chaque modèle est étendu aux mailles voisines les plus
proches. Pour chaque SBIM , une liste des SBIM avec une contribution non nulle
est créée, illustrée par la Figure 41 (d).

- Calcul de dose pour chaque position d’intérêt :

Supposons que est dans le SBIM , la dose déposée en est calculée de la manière
suivante :

Le temps de calcul est exprimé par :

Où est le temps de la propagation des modèles dans les mailles, est le nombre de point d’intérêt et
est le temps moyen pour évaluer la dose déposée en un point donné. Le temps de calcul est quasi-
proportionnel au nombre de points d’intérêt. Il peut être très court pour un nombre de points réduit, ce
qui est très intéressant comme résultat.

Figure 41 Illustration de l’algorithme de calcul de dose

59
 Chapitre IV Contribution Doséclair

Dans ce chapitre, nous expliquons les principaux apports de cette thèse dans la méthode de calcul de
dose. Le calcul de dose se fait par une projection de la maille du fantôme hétérogène dans son modèle
homogène associé. La projection est réexprimée ici par un tenseur pour chaque SBIM. La conception
du tenseur est présentée dans la section 1.

Le processus de propagation du faisceau est caractérisé par la propagation de ces tenseurs de

projection de maille à maille. La propagation à travers les hétérogénéités axiales et latérales du
fantôme constitue la principale difficulté du sujet. Nous automatisons ici le processus en le
réexprimant lui aussi par des tenseurs, ce qui permet au final une définition mathématique unique et
générique du processus de propagation. Nous allons détailler la propagation axiale en section 2, la
propagation latérale en section 3. Enfin, à partir de ces différents éléments, nous présentons quelques
résultats élaborés de ce travail de reformulation.

1 Les projections
Pour chaque SBIM, on associe à sa maille support une projection vers le modèle pré-calculé du même
matériau. Chaque projection est représentée par un tenseur d’ordre 2, une fois variant et une fois
contravariant (il se représente comme une matrice 2D). Ce tenseur agit donc sur un tenseur d’ordre 1
(un vecteur) représentant les coordonnées d’un point du fantôme patient pour donner un autre tenseur
d’ordre 1 représentant les coordonnées de la projection de ce point dans le fantôme homogène. Nous
allons maintenant présenter les repères de coordonnées mis en correspondance par les tenseurs de
projection. Nous regarderons ensuite la forme que prennent ces tenseurs de projection, puis leur calcul
sur la maille support du SBIM.

1.1.Les repères de coordonnées

Afin de mettre en place les tenseurs de projection, il est nécessaire d’associer un repère de
coordonnées à chacune des trois parties mises en correspondance par ces tenseurs, à savoir le fantôme
hétérogène, le faisceau et les fantômes homogènes utilisés dans les modèles pré-calculés.

Dans notre étude, le fantôme hétérogène utilisé est composé de parallélépipède rectangle. Nous
définissons donc le repère comme le repère orthonormé direct dont l’origine est le
centre du fantôme et les axes correspondent aux directions des plans supports des parallélépipèdes
rectangles. Cette définition va nous permettre de simplifier un certain nombre d’écritures et de calculs
(constantes nulles, valeurs d’intérêts correspondant aux projections sur ces axes) sans pour autant
réduire la capacité de généraliser ces calculs sur des modèles plus complexe. Le repère cartésien
est considéré comme le repère de référence. C’est-à-dire que dans la suite de ce mémoire, s’il n’y a pas
de précision spéciale, les vecteurs sont exprimés dans ce repère.

Le repère orthonormé direct , que l’on appelle , est associé à un faisceau

élémentaire. L’origine est le centre (au moins symbolique) de l’origine du faisceau, c’est-à-dire est
le centre de la section de sortie de l’accélérateur. L’axe correspond à l’axe principal du
déplacement des photons primaires, c’est-à-dire l’axe de propagation du faisceau. Les faisceaux
utilisés ayant une section carrée, les axes et sont orientés parallèlement aux bords de la section
de ces faisceaux. Le repère est considéré comme un repère relatif. Dans le processus de projection
que nous souhaitons modéliser, il permet de passer d’un point absolu du fantôme patient à un point
positionné géométriquement par rapport au faisceau, mais sans avoir pris en compte à ce stade la
nature des matériaux et les hétérogénéités.

60
Le repère orthonormé direct , que l’on appelle associé au modèle homogène est
quasiment le même du repère en considérant qu’il n’y a pas de vide entre la sortie (symbolique) de
l’accélérateur et le milieu homogène. L’origine de ce repère est donc à la fois le centre-origine du
faisceau et le centre de la section d’entrée du faisceau dans le matériau homogène. Le milieu étant
anisotrope puisque homogène, le choix de définir le repère en utilisant les axes liés au faisceau (selon
la définition de ) a pour seul objectif et conséquence de simplifier les écritures sans perte de la
capacité de généralisation. Dans le processus de projection que nous souhaitons modéliser, le
« passage » de à permet de se contenter de prendre en compte les hétérogénéités sans avoir à
gérer la géométrie du fantôme patient (qui a déjà été fait dans le passage de à ).

1.2.Le point origine du faisceau RCMI exprimé en coordonnée

sphérique
Le faisceau RCMI est représenté par des faisceaux élémentaires mis côte à côte selon deux directions,
illustré par la Figure 42. Le centre origine du faisceau RCMI est et les faisceaux élémentaires sont
indexés par . Les indices de (resp. ) sont compris entre 1 et (resp. ). Classiquement en
clinique, = =10.

Figure 42 Faisceau RCMI

Le centre-origine du faisceau élémentaire est alors défini par :

,
où est la longueur de côté du faisceau élémentaire (par défaut, la valeur classique dans la littérature
est 1cm).

Figure 43 Coordonnée sphérique

Dans le cadre de notre étude, nous considérons que est placé sur une sphère de rayon fixe, centrée
sur le point centre du fantôme patient et enveloppant ce fantôme hétérogène, illustré sur la Figure
43. peut alors être caractérisé par ses coordonnées sphériques : les angles (longitude) et
(colatitude). Soit le rayon de la sphère (choisi donc de sorte que le fantôme patient soit totalement à
l’intérieur de la sphère (avec même une couronne de sécurité) pour éviter que le faisceau

61
« n’apparaisse » spontanément à l’intérieur du patient), les coordonnées de exprimées dans

le repère sont alors :

Donc le centre du faisceau élémentaire indicé par est défini :

1.3. Formulation d’une projection

Figure 44 Schéma de projection

La Figure 44 présente un cas simple de faisceau et les deux repères d’intérêt pour nous : le repère
fantôme et le repère faisceau . Tout point dans le repère fantôme peut avoir ses coordonnées
exprimées dans le repère faisceau grâce à une fonction de changement de repère sous forme d’un
tenseur d’ordre 2 (matrice) de dimension 4×4:

Soit un point ayant les coordonnées dans le repère fantôme , ses coordonnées dans le
repère faisceau peuvent être déduites par la relation tensorielle : . Cet opérateur
s’appelle un produit tensoriel contracté 1 fois, qui est équivalent à un produit matriciel dans ce cas
présent.

L’orientation du faisceau est caractérisée par le triplet orthonormé de vecteurs . Les trois
premières colonnes du tenseur peuvent être déduites des coordonnées de ces vecteurs exprimés dans
selon :

62
Mais le repère étant orthonormé, les neuf valeurs ne sont pas indépendantes. En fait trois paramètres
suffisent à caractériser l’ensemble des triplets possibles. Il suffit de décrire une succession de rotations
appliquées à un triplet initial (celui de ). Il existe plusieurs façons de décrire ces rotations en 3D.
Les angles d'Euler sont les plus souvent utilisés et nous en avons pris une variante (nommée Tait-
Bryan z-y-x). Ce choix facilite davantage l’optimisation de l’orientation du faisceau que nous allons
expliquer dans la troisième partie de ce mémoire.

Dans la variante choisie, la succession de rotations en défini par :

- une rotation par autour de l'axe ,
- une rotation par autour du nouvel axe (parfois aussi noté )
- une rotation par autour du nouvel axe (parfois aussi noté ).

Les coordonnées du triplet sont alors directement dépendantes de ces trois angles, et pour
la partie du tenseur qui nous intéresse nous obtenons :

Les angles sont bornés :

Pour finir le changement de repère, il reste à prendre en compte le point d’origine du nouveau repère,
ce qui représente une simple translation à partir du repère de référence. Cette translation intervient
dans la définition de la dernière colonne du tenseur et dans l’écriture RCMI directe on obtient :

Grâce à trois angles d’Euler et le centre-origine du faisceau, nous savons passer du repère au
repère grâce à ce tenseur de projection .

Si, de la même façon, on considère un tenseur d’ordre 2 pour passer du repère au repère , en
considérant la particularité de ce dernier par rapport à et le parcours du faisceau, il suffit d’ajuster
par un offset la composante en ce qui correspond à un tenseur de la forme :

Le coefficient traduit les hétérogénéités du milieu, y compris la zone de « vide » entre l’origine de
faisceau et le corps du patient. On peut donc également considérer directement le tenseur de passage
de à , noté . Alors :

63
On note . Nous allons voir plus en détail comment calculer cette composante
par les règles de propagation axiale.

2 Propagation axiale
2.1 Tenseur de la maille support
Chaque SBIM est associé à un tenseur de projection. Ce tenseur permet de projeter tous les points de
la maille support du SBIM dans son modèle homogène. Ce tenseur de projection est propagé
axialement dans les mailles support du fantôme. Dit autrement, pour un SBIM donné, son tenseur de
projection peut être entièrement calculé à partir du tenseur de projection du SBIM qui le
« précède géométriquement ». Comme la détermination des tenseurs à l’entrée du faisceau dans le
fantôme – donc des « tenseurs initiaux » – peut également être établi simplement, nous pouvons donc
évaluer par un calcul de proche en proche l’ensemble des tenseurs. L’objet de cette section est
justement d’établir les « règles » de propagation permettant la systématisation de cette action par un
algorithme pour calculer ces tenseurs.

Pendant la propagation, l’orientation et l’origine du faisceau sont fixes, donc, la « partie »

correspondant à la projection de à (i.e. ne change pas. En revanche, nous allons nous
focaliser sur l’ajustement de la composante qui correspond au passage de à en gérant
justement les hétérogénéités.

De façon générique, pour un point , on chercher à connaitre sa projection dans le repère pour
évaluer la dose qui est pré-calculée dans le modèle en milieu homogène associé au SBIM considéré.
Sa coordonnée selon l’axe , varie de proche en proche en fonction du parcours du faisceau
dans le fantôme hétérogène. L’interprétation à lui donner est qu’elle nous indique à quelle profondeur
dans le modèle homogène nous allons chercher la dose. Cette profondeur s’appelle la profondeur
radiologique. Par définition, pour un point d’un fantôme hétérogène, sa profondeur radiologique,
relativement à un matériau donné est la profondeur où la fluence aurait été autant atténuée dans un
fantôme homogène de ce matériau.

Pour rappel, le principe sous-jacent exploité par la méthode Doséclair, est de considérer que si
localement le milieu et le flux sont homogènes (sur un volume de l’ordre de grandeur des noyaux de
diffusion), même si le fantôme global est hétérogène, la dose localement peut être déduite de la dose
en milieu homogène en prenant la partie soumise à la même fluence. C’est ce principe que traduit
l’utilisation de la profondeur radiologique. Rappelons aussi que lorsque le volume n’est pas assez
grand localement, cela se traduit par la co-influence de plusieurs SBIM qui sont alors combinés par
une somme pondérée selon des mécanismes déjà décrits dans le précédent chapitre.

Nous introduisons maintenant plusieurs concepts et notations pour pouvoir organiser les calculs :
 Afin de simplifier les écritures, nous allons considérer que par rapport à un SBIM considéré
noté , le SBIM noté correspond au SBIM qui « le précède géométriquement » et
que le SBIM noté correspond au SBIM qui « le suit géométriquement »4. Cette

4
Cette notation ne gère pas les découpages ou regroupements de SBIM. On peut toujours considérer que si
plusieurs SBIM sont initiés à partir du SBIM (voire règles de découpage d’un faisceau en SBIM [Deuxième

64
 notation sera utilisée en exposant des fonctions pour les associer aux bons SBIM. Ainsi si on
reprend la Figure 32 et que le SBIM est le SBIM 3 (poumon en bas), alors le SBIM
est le SBIM 1 et le SBIM est le SBIM 5.
 Nous notons le barycentre de la section d’entrée du SBIM dans sa maille support et
le barycentre de la section de sortie du SBIM de sa maille support. Par définition des
SBIM qui précèdent et suivent géométriquement, on a et .5
 La notion de profondeur radiologique d’un SBIM , noté , qui est la profondeur
radiologique du barycentre de la section d’entrée du SBIM dans le modèle homogène du
milieu de sa maille support : .
 La notion de profondeur de SBIM qui correspond à la différence de profondeur
radiologique entre le barycentre de la section d’entrée du SBIM (i.e. la profondeur
radiologique du SBIM) et la profondeur radiologique du barycentre de la section de sortie du
SBIM : . Comme la maille support est homogène, cette
différente de profondeur est strictement géométrique et peut aussi être évaluée dans si
nécessaire par la simple différence d’altitude .

Prenons l’exemple de la Figure 98 de l’annexe A, où il n’y a qu’un SBIM par maille et un découpage
très simple. Si on considère en référence pour notre calcul le second SBIM, alors le SBIM d’entrée
est et le SBIM de sorti est . Si on considère que chaque maille est un milieu différent, il
faut tenir compte de ces changements de milieu dans le calcul de la projection. La profondeur
radiologique du SBIM joue un rôle très important.

Si on regarde la projection du SBIM , elle doit projeter un point du fantôme hétérogène vers un
point du modèle homogène correspondant où la fluence a subi la même atténuation, de façon à assurer
la continuité de la fluence aux interfaces dans le fantôme hétérogène.

Considérons ce qui se passe à l’interface qui sépare les SBIM et . La fluence à l’entrée du
fantôme est . Rappelons-nous que la fluence primaire est atténuée exponentiellement avec la
distance parcourue dans le milieu. Si le milieu possède un coefficient d’atténuation linéique , la
fluence à la profondeur est donnée par : .

En tenant compte de cette notation et en considérant que est le coefficient d’atténuation linéique
du milieu de la maille support au SBIM , alors la fluence au barycentre de la section de sortie du
SBIM vaut :

La fluence au barycentre de la section d’entrée du SBIM est:

Comme il s’agit du même point, grâce à l’hypothèse de la continuité de la fluence aux interfaces :

partie, Chapitre III, Section 3]), alors ils sont indicés en , , etc. De la même façon on pourrait
intégrer pour un SBIM issus de regroupement des et également. D’un point de vue
informatique, dans le simulateur, ces choses sont gérées mais pour la démonstration mathématique afin d’établir
les formules, cette complexification des notations et des situations est inutile.
5
Pour faire écho à la note précédente, on pourrait en cas de découpage avoir plusieurs sections de sortie et
utiliser alors des notations comme , ,…

65
On en déduit facilement que :

La profondeur radiologique du SBIM « successeur » peut donc être calculée à partir du SBIM courant
selon :

De ce fait, nous pouvons en déduire la formule récursive de la composante du tenseur de

projection du SBIM successeur suivant :

+ +

est la distance parcourue depuis le point d’origine du faisceau selon la direction pour
atteindre la profondeur (géométrique) du barycentre d’entrée du SBIM . Cette formule nous
permet de pouvoir propager le tenseur de projection du SBIM au SBIM .

Une validation de cohérence simple à évaluer consiste à considérer le cas trivial : = . Par
conséquent, nous avons égalité entre et , qui implique que le tenseur de projection du
SBIM est propagé tel quel au SBIM suivant sans aucune modification. Ceci est logique puisqu’il
n’y a pas d’hétérogénéité et que si les deux mailles étaient réunies en une seule maille support, les
deux SBIM seraient également réunis en un seul, avec un tenseur unique valide sur toute la grande
maille support.

Le détail de ces calculs est présenté en Annexe A.

2.2 Simulation et résultat

Toutes les évaluations présentées dans ce chapitre se baseront sur un même fantôme muni d’un
maillage régulier de 45 cubes (3×3×5 = 45) de 4 cm de côté, représenté sur la Figure 45. Les doses
sont calculées dans des voxels de 2 mm de côté. Le fantôme test comptera donc 60×60×100 =
45×(4×5)3 = 360000 voxels.

Le faisceau utilisé dans ces simulations est mono énergétique de 5 MeV, il a une fluence homogène
sur toute sa section de forme carrée, dont la surface est de 1×1cm2.

La qualité des résultats de notre méthode est évaluée par comparaison aux résultats obtenus par le code
Monte Carlo Penelope (Salvat, et al., 2006).

Ces comparaisons se feront par l’affichage de la différence de dose en chaque point, et/ou par
l’affichage de l’index delta (Blanpain, et al., 2009) outil équivalent au classique index gamma (Low, et
al., 1998), mais plus rapide à calculer que ce dernier et plus robuste à certains artéfacts de calculs
surtout dans les zones de très fort gradient. Les tolérances en position et en dose retenues pour l’index
delta sont respectivement de 2 mm et de 2%. Le détail de calcul de l’index delta est présenté en
Annexe F.

Nous étudions d'abord le cas d'un faisceau envoyé en plein cœur du fantôme avec 5 couches
composées de différents matériaux, illustré par la Figure 45. Sur ce test, on s’intéresse particulièrement

66
à la dose déposée sur l’axe du faisceau pour vérifier la validité des tenseurs de projection dans les
mailles support.

Sur la Figure 46, la courbe de dose est constituée d’une succession de plateaux décroissants,
correspondant aux différentes mailles. Pour bien illustrer les principes de pondération, nous avons
établi deux types de lois.

 Tout d’abord une loi binaire où le SBIM à une influence de 1 dans sa maille support et 0
ailleurs. On voit que la distribution de dose obtenue au cœur de ces mailles est très précise. Dit
autrement, les mailles supports sont projetées aux bonnes profondeurs dans leurs modèles
homogènes. Par contre, les résultats sont très imprécis près des interfaces, les modifications
lentes des flux électroniques, d’un matériau à l’autre, n’étant pas prises en compte avec cette
loi de pondération.
 Ensuite les fonctions de pondération adaptées que nous avons détaillées dans la section 5 du
chapitre précédent. Grâce à elles, les variations de dose autour des interfaces sont précisément
retrouvées.

Cette précision est confirmée par l’évaluation réalisée via l’index delta et affichée sur la Figure 47. On
voit que l’erreur est très acceptable, puisque pour tout point, delta est compris entre −1 et 1 (et même
entre -0.6 et 0.7).

Figure 45 Un faisceau traversant les cinq couches du fantôme test

Figure 46 La dose déposée sur l’axe du faisceau, Doséclair vs Penelope

67
 Figure 47 Évaluation de la dose sur l’axe du faisceau via l’index delta

3 Propagation latérale
3.1 Tenseur de la maille voisine
Nous avons créé les tenseurs de projection pour les mailles support du fantôme. Ensuite, ces tenseurs
sont utilisés pour projeter la maille support dans leur modèle homogène. Du fait de la diffusion latérale
par les électrons (et ainsi, dans une moindre mesure, par les photons secondaires), la dose peut se
déposer à plusieurs centimètres de chaque côté du faisceau primaire. Dans cette section, nous allons
voir comment créer les tenseurs de projection dans les mailles voisines.

Pour une maille voisine donnée, on cherche à calculer le tenseur de projection qui permet de projeter
les points de cette maille voisine dans le modèle homogène de même matériau que la maille support.
Pour cela, nous allons étendre les tenseurs de maille support à leur mailles voisines, ce que l’on
appelle la propagation latérale du tenseur. Cette façon de faire a été étudié et juger la meilleure parmi
plusieurs approche dans la thèse de (Blanpain, et al., 2009).

Figure 48 profils de dose avec un faisceau de photons de 5 MeV, de section carrée 1×1 cm2

Sur la Figure 48, le faisceau longe une interface (représentée par le trait bleu) entre deux mailles. On
voit deux profils de dose en fonction de l’écart à l’axe latéral du faisceau.

Dans le premier cas, les deux mailles sont composées d’eau. Le profil est présenté par la courbe en
pointillé noir. Or la maille voisine est composée du même matériau que la maille support, la diffusion
latérale se fait naturellement, exprimée par une continuité de la courbe dans la maille voisine. Ceci
s’est traduit par la propagation directe du tenseur de la maille support à sa maille voisine.

68
Dans le deuxième cas, la maille à droite est composée d’os. Il apparait évident qu’un changement de
matériau dans la maille voisine peut avoir d’importantes conséquences sur la décroissance latérale de
la dose. On peut voir que la dose décroit plus vite dans l’os (courbe pleine rouge) que dans l’eau. Ceci
est dû au fait que la densité électronique de l’os est supérieure à celle de l’eau et est également
cohérent avec une valeur du coefficient linéique plus importante. Cette hétérogénéité latérale du
fantôme impacte sur la distance latérale projetée dans le modèle homogène de la maille support.

Figure 49 Propagation latérale avec l’hétérogénéité

Ce fantôme hétérogène est représenté par la Figure 49. Le tenseur de projection du SBIM de la
maille d’eau est propagé dans la maille voisine d’os par l’interface encadrée en pointillé rouge. La
dose déposée au point dans la maille d’os est égale à la dose déposée au point dans le modèle
homogène d’eau. On note la projection orthogonale du point sur l’interface. Alors le tenseur
de projection de cette maille voisine, qui permet de projeter les points de la maille d’os dans le modèle
homogène d’eau, doit modifier (en fait ici augmenter) la vraie distance latérale géométrique
selon un coefficient de scaling pour que la dose soit égale à la bonne dose à la nouvelle distance
. Cette nouvelle distance est appelé la distance latérale projetée. C’est ce scaling partiel qui
va permettre de projeter le point de la maille d’os au point dans le modèle homogène d’eau.
Cette projection est intégrée dans la définition du nouveau tenseur de la maille d’os.

Pour le choix du coefficient de scaling, (Blanpain, et al., 2009) a testé avec le ratio de coefficient
d’atténuation linéique, mais cela ne fournit pas de résultats satisfaisants. Ces ratios donnent en effet
des scalings trop élevés, ou trop faibles, selon les cas. (Blanpain, et al., 2009) a donc choisi d’utiliser
des coefficients de scaling obtenus expérimentalement, en utilisant, pour chacun des milieux
homogènes simulés, la distance entre le bord du faisceau et le point où la dose est de 20% de la dose
sur l’axe du faisceau. On obtient ensuite les coefficients en faisant des rapports entre ces valeurs.

Figure 50 Zone de validité pour le tenseur latéral

Un exemple plus complexe est présenté par la Figure 50, le faisceau traverse obliquement le fantôme.
Dans la maille d’eau en bas à gauche, un SBIM en pointillé rouge est créé à l’entrée du fantôme. Ce
SBIM est propagé dans la maille poumon en axial et latéral par l’interface 1. Ici, on ne s’intéresse qu’à
la propagation latérale du tenseur dans la maille poumon. En pratique, la maille poumon est divisée en
deux volumes distincts qui sont séparés par la surface en pointillé noir. À chacun de ces deux volumes,
un tenseur de projection latérale différent est associé. Précisément, pour le volume en haut (resp. en

69
bas) de la maille, le scaling est appliqué sur la distance latérale par rapport à l’interface 2 (resp.
interface 1). La création des interfaces de scaling candidates est matérialisée dans la propagation
latérale du tenseur.

Figure 51 Propagation latérale multiple

Enfin, une propagation latérale multiple est aussi envisageable à l’aide du principe de propagation de
maille à maille. Sur l’exemple de la Figure 51, le SBIM est propagé en latéral d’abord dans la
maille d’os, ensuite dans la maille poumon. Le tenseur du SBIM est propagé dans la maille d’os, et
ce nouveau tenseur est ensuite propagé dans la maille du poumon. Chaque propagation est réalisée par
une addition tensorielle.

En pratique, la propagation du tenseur de projection est bien plus complexe, en particulier dans
l’espace à trois dimensions. Les approximations et la réalisation de cette propagation sont détaillées en
Annexe B.

3.2 Simulation et résultat

Dans cette section, nous allons vérifier la propagation latérale du tenseur de projection. Nous avons
repris les tests présentés dans les travaux de (Blanpain, 2009) et recalculé la dose déposée au travers
de la formulation tensorielle.

Figure 52 Un faisceau traversant le fantôme en longeant au plus près des interfaces

70
Figure 53 Simulation faisceau longeant les interfaces latérales. Le plan de coup de la distribution de dose calculée par
Penelope (a), par Doséclair (b). Évaluation de l’index delta de la dose Doséclair comparée à Penelope (c).

D’abord, nous étudions le cas d'un faisceau horizontal traversant les cinq mailles centrales mais en
étant accolé aux interfaces latérales, illustré en Figure 52. La distribution de dose obtenue est identique
à celle de (Blanpain, 2009). Nous pouvons voir sur la Figure 53(a) et (b), la diffusion latérale à
proximité des deux hétérogénéités est maîtrisée, la dose calculée par Doséclair ne diffère de celle de
Pénélope que sur les bords du faisceau où la dose n’est plus significative par rapport à celle déposée à
l’axe du faisceau. Ce résultat est confirmé par l’évaluation de l’index delta (Figure 53(c)), les écarts de
dose se situent au voisinage des interfaces, mais la valeur de l’index delta de ces zones est comprise
entre -1 et 1.

Figure 54 Un faisceau traversant obliquement le fantôme

71
Figure 55 Simulation faisceau oblique. Le plan de coup de la distribution de dose calculée par Penelope (a), par
Doséclair (b). Évaluation de l’index delta de la dose Doséclair comparée à Penelope (c). Les zones en dehors des
limites tolérées (d).

Ensuite, nous étudions le cas d'un faisceau oblique traversant le fantôme, illustré en Figure 54. La
distribution de dose obtenue est identique à que celle de (Blanpain, 2009). Nous pouvons voir sur la
Figure 55 (a) et (b) que la dose déposée à l’entrée du fantôme est très proche de la dose Penelope
malgré son entrée oblique. En revanche, des erreurs sont non négligeables pour certaines positions
situées à proximité des interfaces eau/os et eau/poumon, comme montré par l’évaluation de l’index
delta (Figure 55 (c)). Les voxels où l’erreur dépasse la tolérance sont colorées en rouge, illustré par la
Figure 55 (d). Sur ce plan de coupe (6000 voxels), 21 voxels sont en dehors du critère d’acceptabilité,
qui représente un taux de 0.35%. 154 voxels le sont pour le fantôme entier, qui représente un taux de
0.042%.

4 Discussion
Les simulations présentées dans cette section utilisent un faisceau mono-énergétique de 5Mev. Notons
que les spectres en énergie du faisceau de photon sortant des accélérateurs linéaires de différents
modèles sont compris entre 500keV et 18 MeV. Dans un faisceau poly-énergétique, les photons de
basse énergie sont éliminés plus rapidement que ceux de forte énergie. L’énergie moyenne des photons
du faisceau augmente avec la profondeur. C’est ce qu’on appelle le durcissement du faisceau. Ce
durcissement du faisceau a des conséquences sur la dose déposée. Un faisceau durci transmet une
énergie en moyenne plus élevée aux électrons qui auront des trajectoires plus longues, modifiant ainsi
les distributions de dose (Metcalfe, et al., 1990) (Bjarngard, et al.). La méthode Doséclair gère ce
phénomène de façon transparente en utilisant directement le dépôt de dose du même faisceau poly-
énergétique dans les modèles homogènes. Les propriétés particulières de la fonction exponentielle et le
fait de travailler aux interfaces avec les profondeurs équivalentes entre les milieux assurent une bonne
continuité également du spectre implicite. Quant à la propagation, le coefficient d’atténuation et le

72
scaling latéral ne sont pas affectés par la variation de l’énergie moyenne des photons car pour la plage
d’énergie entre 500keV et 18 MeV, les atténuations des rayons sont dominées par l’effet Compton
(énergies supérieures à 110 keV). Pour les interactions de type Compton, la probabilité d’interaction
( ) dépend essentiellement de la densité électronique de la matière (nature du milieu) et peu de
l’énergie des photons. Donc, l’atténuation des rayons poly-énergétiques par effet Compton dépend
principalement de la nature du milieu, qui est indépendante de l’énergie du faisceau. Sachant que le
principe du modèle de propagation des SBIMs repose sur une mise à échelle selon les natures des
milieux, le modèle de la propagation n’est pas affecté par la présence des rayons poly-énergétiques.

Avec la reformulation tensorielle, nous sommes capables de reproduire la même distribution de dose
que (Blanpain, 2009). Bien que n’apportant rien de nouveau sur la connaissance scientifique même du
modèle, cette reformulation tensorielle a trois grands avantages.

Tout d’abord, elle permet de couvrir l’ensemble des situations géométriques en un calcul et une
formulation finale unique. Par rapport aux formulations matricielles plus directes de (Blanpain, 2009)
qui donnait une formule pour chaque situation géométrique, ne se pose plus ici la question de
l’exhaustivité des situations envisagées, que ce soit pour les calculs ou pour l’implémentation
informatique. Cette dernière est d’ailleurs considérablement compactée.

Cette reformulation ouvre aussi la voie au calcul du gradient de dose et donc à l’optimisation qui font
l’objet de la partie suivante. Cette perspective à elle seule d’ailleurs justifie ce travail de reformulation.
Et de la même manière, en définissant le gradient et la propagation du gradient par des tenseurs, on
obtient une démonstration et une formulation finale plus compliquées mathématiquement, mais unique
(on ne refait pas les calculs pour chaque cas particulier) et beaucoup plus simple à implémenter
informatiquement.

Enfin cette reformulation tensorielle ouvre aussi la voie à des fonctions de projections plus complexes.
En effet en « agrandissant » le tenseur de position (équivalent à dans le cas présenté) avec des
grandeurs comme , , , , et , la même écriture tensorielle permet alors de faire des
projections selon un polynôme d’ordre deux. Une telle projection est l’une des pistes envisagées pour
identifier les formules de projection et de propagation lorsque les mailles sont définies par des
quadriques (cas que l’on trouve dans certain outils de Monte-Carlo et qui permet d’utiliser beaucoup
moins de mailles pour définir les organes). Les heuristiques restent encore à être définies bien sûr,
mais une fois établies, toutes les démonstrations faites ici ainsi que les implémentations logicielles
sont à légèrement adapter mais restent valides.

73
 Troisième partie

Optimisation du plan de traitement

74
Chapitre I Méthodes d’optimisation de traitement : un état de l’art ............................................... 77
1 Les critères d’optimisation ............................................................................................................ 78
2 Formulation mathématique du problème ....................................................................................... 78
3 Optimisation des intensités seules ................................................................................................. 80
3.1 Modèle linéaire ...................................................................................................................... 80
3.2 Modèle non linéaire ............................................................................................................... 80
4 Optimisation des incidences .......................................................................................................... 81
4.1 Méthodes exhaustives............................................................................................................ 82
4.2 Méthodes itératives déterministes ......................................................................................... 82
4.3 Méthodes itératives stochastiques ......................................................................................... 83
Chapitre II Optimisation des plans de traitement en RCMI............................................................. 85
1 Modèle optimisation et formulation mathématique....................................................................... 85
1.1 Objectifs d’optimisation ........................................................................................................ 85
1.2 La modulation d’intensité du faisceau ................................................................................... 86
1.3 Balistique du traitement......................................................................................................... 87
2 Gradient de dose ............................................................................................................................ 87
2.1 Composition du gradient de dose .......................................................................................... 87
2.2 Le gradient du tenseur de projection ..................................................................................... 90
3 Descente de gradient ..................................................................................................................... 90
3.1 Algorithme générale .............................................................................................................. 90
Chapitre III Implémentation logicielle et validation des résultats .................................................... 92
1 Implémentation logicielle .............................................................................................................. 92
2 Validation des paramètres individuels........................................................................................... 93
2.1 Conditions de simulation ....................................................................................................... 93
2.2 Validation du positionnement du point d’origine du faisceau ............................................... 94
2.3 Validation de l’orientation du faisceau.................................................................................. 96
2.3.1 Paramètre ................................................................................................................... 97
2.3.2 Paramètre ................................................................................................................. 100
2.3.3 Paramètre ................................................................................................................. 103
3 Fantôme homogène d’eau avec cible de forme concave ............................................................. 105
3.1 Conditions de simulation ..................................................................................................... 105
3.2 Cyberknife - faisceaux élémentaires indépendants ............................................................. 106
3.2.1 Intensité uniforme (49 mins) ....................................................................................... 106
3.2.2 Intensité variables (24.5 mins) .................................................................................... 110
3.3 Faisceaux RCMI classique .................................................................................................. 114
3.3.1 Plan angles fixes avec modulation (3.7 mins) ............................................................. 114
3.3.2 Plan isocentre (23 mins) .............................................................................................. 117
3.3.3 Plan 3D libre (76 mins) ............................................................................................... 118
4 Fantôme hétérogène Penelope Poumon-Os ................................................................................. 119

75
4.1 Conditions de simulation ..................................................................................................... 119
4.2 Échantillonnage des points de contrôles.............................................................................. 120
4.3 Étude de cas (82 mins) ........................................................................................................ 120

76
 Chapitre I Méthodes d’optimisation de
traitement : un état de l’art
Pour les traitements conventionnels (non modulés), la planification est dite directe et il revient à
l’utilisateur d’optimiser la géométrie d’irradiation (nombre et orientations des faisceaux, limites du
champ d’irradiation, avec ou sans filtre en coin, pondération des faisceaux, etc.) qui permet de
satisfaire les contraintes dosimétriques fixées, illustré en Figure 56 (a). Dans les cas complexes,
l’utilisateur peut être amené à modifier plusieurs fois les paramètres d’irradiation jusqu’à atteindre un
plan de traitement jugé satisfaisant. La procédure «essai et erreur» est réalisée par le planificateur
(dosimétriste).

Dès la mise en place de la RCMI dans les années 90, cette méthode d’essais successifs se révéla
impossible à faire manuellement car le nombre de possibilités de modulation est infini. Il faut donc
faire appel à des algorithmes automatiques de recherche de solutions que l’on appelle l’optimisation
par planification inverse. La planification inverse consiste à déterminer, à partir des objectifs cliniques
définis par le médecin oncologue (dose à atteindre pour la tumeur et dose à ne pas dépasser pour les
organes à risques), les faisceaux modulés à utiliser, illustré en Figure 56 (b). Cette planification diffère
radicalement de celle de conventionnelle car la procédure «essai et erreur» cette fois ci est réalisée par
le calculateur (logiciel).

(a) (b)

Figure 56 Optimisation conventionnelle (a) et optimisation par planification inverse (b) (Lefkopoulos, et al., 1999)

Ce calculateur appelé TPS est basé sur un algorithme de calcul de dose associé à un algorithme
d’optimisation (Wu, et al., 2000) (Brahme, 1988) (Spirou, et al., 1998). À l’heure actuelle6,
l’optimisation dans le TPS ne concerne que la modulation, l’utilisateur doit déterminer au préalable le
nombre et l’orientation des faisceaux. C’est-à-dire que les intensités de chaque faisceau choisi sont
calculées automatiquement pour obtenir la distribution de dose souhaitée. Pour cela, ces algorithmes
traduisent mathématiquement le problème à l’aide d’une fonction objectif exprimant généralement la
différence entre la distribution de dose prescrite et la distribution de dose calculée. Leur but est donc
de minimiser cette fonction objectif.

Le problème de planification d’inverse est un problème complexe et non intuitif, auquel de nombreux
auteurs ont essayé d’apporter des solutions grâce à l’utilisation de méthodes d’optimisation très
variées.

6
Des systèmes permettant l’optimisation simultanée de la géométrie et de la modulation commencent à être
disponibles sur les systèmes commercialisés.

77
1 Les critères d’optimisation
La première grande difficulté du problème d’optimisation réside dans la traduction des objectifs
cliniques sous forme de critères mathématiques afin que la fonction objectif représente au mieux le
problème réel.

Deux types de critères (ou contraintes) peuvent être utilisés : les critères physiques, qui sont basés sur
des contraintes de dose particulières et les critères biologiques, qui sont basés sur les effets biologiques
entraînés par l’irradiation.

Les critères physiques se présentent sous la forme de doses minimales à atteindre et de doses
maximales à ne pas dépasser pour les volumes cibles, et de doses maximales à respecter pour les
OARs7. De tels critères s’appliquent donc principalement aux volumes cibles et aux OARs présentant
une réaction à seuil (organe à morbidité sévère comme la moelle épinière par exemple). Pour les
organes qui ont des volumes importants (les parotides ou les poumons par exemple), des critères
concernant la relation dose-volume ont été introduits. Ils peuvent s’énoncer ainsi : « x% du volume ne
doit pas recevoir une dose supérieure à y Gy » (Spirou, et al., 1998).

Bien que ces critères physiques soient utilisés dans la plupart des systèmes de planification inverse
actuels, pour certains auteurs (Brahme, 1999) (Mohan, et al., 1994) (Wang, et al., 1995), ils ne
suffisent pas à obtenir des résultats totalement satisfaisants d’un point de vue clinique. En effet,
l’indication d’une contrainte de type dose-volume ne renseigne pas suffisamment sur les effets
biologiques entraînés par l’irradiation. Il a alors été proposé d’utiliser des critères biologiques destinés
à quantifier la réussite d’un traitement.

Les critères biologiques peuvent être basés sur les notions de probabilité de contrôle tumoral (« Tumor
Control Probability », TCP) et de probabilité de complication des tissus sains (« Normal Tissue
Complication Probability », NTCP), qui prédisent respectivement le taux de survie des cellules
cancéreuses, et l’impact de la distribution de dose sur les tissus sains environnants. La qualité d’un
plan de traitement sera donnée par le couple TCP/NTCP. Idéalement, on vise TCP = 1 et NTCP = 0.

Le concept de la dose équivalente uniforme (Equivalent Uniform Dose, EUD) introduit par
(Niemierko, 1997) est utilisé pour évaluer la qualité d’un plan de traitement à partir de critères
biologiques. L’EUD permet, en traduisant une distribution de dose inhomogène en une distribution de
dose homogène entraînant le même taux de mort cellulaire, de quantifier la distribution de dose reçue
dans une structure donnée en tenant compte de sa sensibilité aux rayonnements. La formule générale
de l’EUD, valable à la fois pour les tissus tumoraux et les tissus sains, est la suivante :

représente le nombre de voxels dans le volume d’intérêt, est la dose dans le ième voxel et est un
paramètre spécifique du volume d’intérêt qui décrit l’effet dose volume.

Ce paramètre peut être déterminé empiriquement à partir des données dose-volume publiées par
(Emami, et al., 1991). L’avantage de l’EUD est de permettre la conversion d’une distribution de dose
en un seul paramètre prenant en compte un Histogramme Dose-Volume (HDV) dans sa totalité
(contrairement aux contraintes dose volume) ainsi que la réponse tissulaire. Son inconvénient principal
réside dans la difficulté d’établir des valeurs fiables pour le paramètre .

2 Formulation mathématique du problème

7
On pourrait avoir plusieurs PTV avec des seuils différents et plusieurs OAR avec des sensibilités différentes.

78
La modélisation du processus de la RCMI est schématisée sur la Figure 57. Nous présentons pour des
raisons de simplification une configuration avec une distribution de dose en 2D et des faisceaux
modulés en 1D. Cette formulation peut être facilement étendue vers une configuration avec une
distribution de dose en 3D et des faisceaux modulés en 2D.

Figure 57 La modélisation du processus d’irradiation en radiothérapie d’intensité modulée (Lefkopoulos, et al., 1999)

L’équation de calcul de dose (pour construire la distribution de dose) s’exprime par le système
d’équations linéaires suivant :

Où :
est un vecteur de dose, représente la dose déposée dans le ème voxel
est une matrice de contribution de dose, représente la contribution en dose du faisceau
élémentaire uniforme dans le voxel
est un vecteur des intensités, représente l’intensité du faisceau élémentaire

Chaque voxel appartient à un volume d’intérêt. Selon ces volumes, la matrice est partitionnée en
sous-matrices , la lettre ‘c’ désigne le volume cible, ‘OAR’ désigne l’organe à risque.
Chaque matrice représente la contribution de dose dans le volume en question par chaque faisceau. En
tenant compte de cette partition, l’objectif et les contraintes idéales du traitement s’expriment par :

Où , et sont les seuils de dose dans le volume cible et l’OAR respectivement.

Pour définir les éléments de la matrice , différents modèles de calcul de dose peuvent être utilisés.

Rappelons que l’objectif de la planification est de trouver, à partir d’une prescription de dose et des
contraintes, les faisceaux modulés à utiliser. Ces faisceaux modulés sont définis par leurs incidences et
leurs intensités. Le processus de planification est généralement composé de deux étapes
séquentiellement :
1. Le problème de la géométrie où on cherche à trouver les incidences des faisceaux
2. Le problème des intensités où on cherche à déterminer les intensités des faisceaux pour une
incidence retenue

79
Une fois que les incidences des faisceaux sont choisies, les intensités des faisceaux sont calculées.

3 Optimisation des intensités seules

Les objectifs cliniques sont traduits par des fonctions objectifs qui dirigent le processus
d’optimisation. Les deux types de critères (physique et biologique) que l’on a détaillés dans la section
1 de ce chapitre sont utilisés dans la définition de ces fonctions.

3.1 Modèle linéaire

Les fonctions objectifs sont apparues en premier sous forme linéaire. (Hodes, 1974) (Langer, 1987)
(Legras, et al., 1982) (Morrill, et al., 1991) (Powlis, et al., 1989) (Rosen, et al., 1991). Elles utilisent
des critères physiques pour atteindre des objectifs du type : minimiser l’écart moyen entre la dose
déposée et la dose prescrite, minimiser la dose maximale/moyenne dans l’OAR. Les contraintes sont
du type : les intensités sont non-négatives, des seuils de dose ou de dose moyenne sur les volumes
d’intérêt.

L’avantage de la méthode linéaire est sa rapidité à trouver une solution. L’algorithme du simplexe est
souvent utilisé pour le problème linéaire. L’idée est de parcourir le polyèdre convexe formé par les
contraintes linéaires du problème en optimisant la fonction objectif. En revanche, la solution trouvée
par cet algorithme se trouve toujours à l’extrémité du modèle. Concrètement, une solution issue de la
méthode du simplexe génère un plan dans lequel les seuils de dose sont atteints dans le volume cible et
les OARs. En générale, ceci représente un plan sous-optimal. En outre, avec des contraintes sévères, le
modèle linéaire classique ne peut pas toujours admettre une solution. Pour surmonter ce problème de
faisabilité, un nouveau modèle linéaire ayant des contraintes élastiques a été proposé par (Holder,
2003)(Holder, 2004)(Holder, 2006). Les contraintes dose-volume ne sont pas incluses dans ce modèle.

3.2 Modèle non linéaire

Toujours avec des critères physiques, le modèle non-linéaire basé sur les moindres carrées pondérés
est le plus utilisé dans les formulations :

Où :
est la dose prescrite à la cible,
et sont les paramètres de priorité des volumes d’intérêt,
et représentent respectivement les nombre de voxels du volume cible et de l’OAR,
l’opérateur .

À ce jour, cette formulation quadratique est reconnue comme un modèle standard et implémentée dans
plusieurs TPS commerciaux.

L’inconvénient de ce modèle est son choix arbitraire des priorités (via et ). Un processus de tests-
erreurs est imposé pour trouver les bons paramètres afin d’optimiser le plan final.

En réalité, deux distributions de dose différentes peuvent conduire à un effet équivalent vis-à-vis du
contrôle tumoral et de la protection du tissu sain. Une fonction objectif basée sur des critères
biologiques permet aussi de qualifier le plan de traitement.

80
Dans (Kallman, et al., 1992) (Lof, 2000) (Mageras, et al., 1993) (Morrill, et al., 1991) (Wang, et al.,
1995) (Webb, 1992), le modèle est fondé sur les critères TCP (tumor control probabilities) et NTCP
(normal-tissue complication probabilities) que nous avons vu dans le précédent chapitre.

En supposant que l’égalité entre deux probabilités du contrôle tumoral implique l’équivalence des
deux distributions de dose, le critère EUD (tumor-applicable equivalent uniform dose) est utilisé pour
définir la fonction objectif (Niemierko, 1997) (Niemierko, 1999) (Thieke, et al., 2002) (Wu, et al.,
2002).

Les fonctions objectifs biologiques fournissent une mesure de la qualité de vie du patient après
traitement par irradiation. Elles nécessitent des paramètres radiobiologiques qui décrivent les réponses
à la dose des différents volumes d’intérêt. Ce type de modèle a besoin d’être correctement vérifié pour
être fiable et applicable en clinique.

Deux familles d’algorithme ont été utilisées pour résoudre les problèmes non-linéaires en
radiothérapie : les méthodes itératives et les méthodes d’inversion. Un panorama de ces techniques est
proposé par (Webb, 1993).

Les méthodes d’inversion représentent une manière directe de détermination de la solution optimale à
partir de la matrice telle que les techniques de décomposition en valeurs singulières (Platoni, 1998)
(Grandjean, 1998) et CLS (constrained least-squares) (Starkschall, 1984) (Crooks, et al., 2002).

Les méthodes itératives sont les plus souvent utilisées. Elles sont subdivisées en deux catégories : les
méthodes déterministes (descente du gradient) et les méthodes stochastiques (recuit simulé et
algorithmes génétiques, dites aussi métaheuristique).

La descente du gradient dont la direction de recherche est bien définie est relativement rapide en temps
de calcul. Cependant, il y a le risque de tomber dans des minima locaux.

Contrairement à la descente du gradient, les méthodes stochastiques intègrent des mécanismes

aléatoires d’exploration de l’espace de recherche. Ces méthodes présentent l’avantage de converger
vers l’optimum global, mais la vitesse de convergence n’est pas maitrisée. Leur utilisation est plus
adaptée en présence de fonctions objectifs ayant beaucoup de minima locaux, comme les fonctions
objectifs biologiques.

La quasi-totalité des systèmes de planification actuels utilisent la technique du gradient ou une de ses
améliorations, car cette méthode représente le meilleur compromis entre l’obtention d’une solution
satisfaisante et un temps de calcul raisonnable. Bien que cette technique ne permette pas de garantir
que la meilleure solution soit atteinte, il a été montré que, dans le cas d’une optimisation des intensités
seules (sans ajouter l’orientation des faisceaux), les solutions obtenues sont proches de la solution
idéale pour les raisons suivantes (Bortfeld, 1998) :

1. pour des fonctions objectifs simples, on peut montrer qu’il n’existe pas de minimum locaux
2. le choix d’un point de départ proche du minimum global permet d’éviter les minimums locaux
3. les minimums locaux ne sont pas très différents des minimums globaux

4 Optimisation des incidences

Le choix des incidences des faisceaux dans la planification inverse est crucial. Un petit nombre de
faisceaux avec des incidences optimales peut donner de meilleurs résultats qu’un plus grand nombre
de faisceaux avec des incidences sous-optimales (Bortfeld, et al., 1993) (Soderstrom, et al., 1995)
(Stein, et al., 1997).

81
L’interdépendance entre les incidences et les intensités des faisceaux font de ce problème
d’optimisation des incidences des faisceaux, un problème complexe dont l’espace des solutions est très
important.

Bien que de nombreux auteurs se soient attachés à la résolution de ce problème en proposant

différentes méthodes, il reste toujours un problème ouvert, puisqu’à ce jour aucun TPS commercialisé
ne permet l’optimisation automatique des incidences des faisceaux en RCMI.

4.1 Méthodes exhaustives

Ces méthodes consistent à tester toutes les combinaisons possibles d’incidences de faisceaux afin de
déterminer la balistique optimale. Plus le nombre de faisceaux utilisé sera grand, plus les
combinaisons seront nombreuses.

Cette approche exhaustive a notamment été utilisée par (Wang, et al., 2004). Elle consistait à explorer
toutes les combinaisons possibles de 3, 5 et 7 faisceaux, parmi un choix de 36, 24, 19 et 18 incidences
de faisceaux equi-espacés, en utilisant un générateur de combinatoires ( est le nombre d’incidence,
est le nombre de faisceaux) pour explorer toutes les possibilités. Pour se donner un ordre d’idée, si
l’on veut examiner toutes les combinaisons possibles existant entre 3 faisceaux pour 36 incidences par
faisceau, il faudra tester combinaisons.

Ensuite, pour chaque balistique testée, les intensités ont été optimisées par un algorithme de gradient.
Pour que le calcul soit réalisable dans un temps raisonnable, une méthode de calcul de dose simplifiée
et une parallélisation du processus est utilisé.

(Potrebko, et al., 2007) ont également utilisé une méthode exhaustive afin de déterminer l’incidence
optimale du premier faisceau d’une balistique équi-espacée, dans le cas des traitements de la prostate
en RCMI. La particularité de cette méthode est qu’elle se base sur un critère géométrique qui concerne
un OAR à la fois. Cette méthode est simple à mettre en œuvre. Cependant, la couverture du volume
cible n’est pas prise en compte.

Les méthodes exhaustives sont très coûteuses en temps de calcul. La plupart des auteurs ont donc eu
recours à des méthodes de recherches non-exhaustives, itératives, visant à réduire le temps
d’optimisation pour être compatible avec une utilisation en routine clinique.

4.2 Méthodes itératives déterministes

(Gaede, et al., 2004) ont proposé une méthode heuristique pour sélectionner séquentiellement des
faisceaux optimums. Cette méthode a montré des améliorations en termes d’uniformité de couverture
du volume cible et en termes de réduction de dose aux OAR avec un nombre de faisceaux inférieur à
celui utilisé par les plans équi-espacés. Cependant, les faisceaux optimums obtenus sont très proches
dans l’espace. Cela pourrait entraîner un surdosage à la peau et aux tissus sains.

(Lee, et al., 2006) ont suggéré d’utiliser la programmation linéaire mixte pour optimiser
simultanément les orientations et le nombre de faisceaux ainsi que les intensités. Les variables binaires
sont utilisées pour la sélection des faisceaux, les variables réelles positives sont utilisées pour définir
les intensités. Les critères de dose et de dose-volume sont redéfinis par des contraintes d’inégalités
linéaires. Cette méthode permet, grâce à la redéfinition de la fonction objectif et des contraintes,
d’augmenter le nombre de critères d’optimisation, et de les diversifier. Cependant, elle reste une
méthode de sélection à partir d’un ensemble des faisceaux candidats coplanaires.

(Craft, 2007) a par ailleurs proposé d’optimiser les incidences par une méthode de gradient dans le
cadre d’une formulation linéaire. L’optimisation s’est faite en local sur des incidences initialisées
aléatoirement. La recherche locale a permis une amélioration de la fonction objectif de 6% sur un
fantôme test représentant un cancer du pancréas. Cependant, cette méthode a besoin d’être prouvée

82
pour des formulations plus complexes – i.e. non linéaires – des critères à optimiser. De plus, les
incidences étudiées restent coplanaires. La méthode oblige un stockage de la dose déposée par chaque
incidence discrète. Pour une configuration en 3D, ceci augmente considérablement le besoin de
l’espace mémoire.

4.3 Méthodes itératives stochastiques

L’algorithme du recuit simulé a très souvent été cité dans la littérature pour optimiser les incidences de
faisceaux dans le cas des traitements par RCMI. L’algorithme nécessite toutefois des temps de calcul
très importants, non compatibles avec une utilisation en routine clinique. De nombreux auteurs ont
essayé d’améliorer la vitesse de convergence vers la solution optimale, par différents moyens.

(Djajaputra, et al., 2003) ont utilisé un recuit simulé accéléré pour optimiser les incidences des
faisceaux. L’accélération est réalisée grâce à une loi de décroissance rapide de la température et un
kernel compact pré-calculé pour chaque angle du faisceau. Le temps de calcul a été fortement réduit
pour trouver le plan optimal sur le cas de la prostate et l’ORL.

(Rowbottom, et al., 2001) ont introduit des critères a priori visant à réduire l’ensemble des faisceaux
candidats. Ces critères sont basés sur la dose déposée à l’OAR (qui a une faible tolérance aux
irradiations) dans le but de ne pas avoir de faisceaux traversant ces derniers. Chaque faisceau est
évalué individuellement à l’aide d’une fonction objectif en tenant compte de ce critère. Les auteurs ont
montré que l’utilisation de 3 faisceaux d’orientations optimisées a permis d’améliorer les résultats
obtenus avec 7 (ou 9) faisceaux équi-espacés, dans le cas d’une tumeur de la glande parotide.

(Pugachev, et al., 2002) (Pugachev, et al., 2001) ont réalisé une élimination au préalable de faisceaux
en se basant sur le score BEVD (Beam’s Eye View Dosimetric). L’utilisation du critère BEVD permet
de s’assurer que tous les faisceaux utilisés couvrent suffisamment le volume cible et que les tolérances
des OARs ne sont pas dépassées. Ce critère semble bien représenter l’objectif clinique mais pourrait
être amélioré en incluant des contraintes dose-volume.

Cependant, ces deux méthodes basées sur des critères de pertinence a priori analysent le score de
chaque faisceau individuellement. L’interdépendance des faisceaux constituant la balistique n’est donc
pas prise en compte.

(Wang, et al., 2004) ont montré que l’utilisation du seul critère BEVD ne conduit pas à la même
sélection des angles que ceux choisis par la technique exhaustive en testant toutes les combinaisons
existantes. Par conséquent, l’utilisation de critères évaluant la qualité de chaque faisceau sans tenir
compte de l’influence des autres faisceaux constituant la balistique, ne semble pas adaptée pour
déterminer la balistique optimale.

Ainsi, les algorithmes génétiques ont également été grandement appliqués dans l’optimisation des
incidences des faisceaux. (Hou, et al., 2003) (Li, et al., 2004) (Schreibmann, et al., 2004) ont montré
que ces algorithmes ont permis d’améliorer les distributions de doses obtenues manuellement par
essais successifs ou en utilisant des faisceaux équi-espacés. Cependant, du fait de la taille très
importante du problème, le processus nécessite un temps de calcul très long malgré les qualités de
cette métaheuristique sur ce point. De plus, la détermination de la solution optimale dépendra donc
toujours très fortement de la complexité du problème et des paramètres de l’algorithme.

Les deux aspects du problème de planification du traitement ont été présentés. Le problème des
intensités est d’avantage abordé depuis la mise en place de la RCMI. Différents algorithmes sont
proposés selon la formulation de la fonction objectif et les contraintes. Le problème de la géométrie
n’est pas réellement résolu puisque toujours fait manuellement en clinique. Ce problème est difficile à
résoudre à cause de sa dépendance avec le problème des intensités. Le couplage des deux problèmes
agrandit considérablement l’espace des solutions. Les méthodes existantes pour résoudre les deux

83
problèmes simultanément effectuent une sélection des faisceaux discrets préalablement définis à l’aide
d’une formulation linéaire.

Dans le prochain chapitre, nous allons présenter une nouvelle méthode permettant d’optimiser les
intensités et les incidences 3D de façon continue en les paramètres. L’optimisation se fait par descente
du gradient, simultanément sur ces deux éléments. Ensuite, nous allons montrer la faisabilité de cette
méthode sur des cas synthétique et les résultats révèlent une amélioration de la distribution de dose par
rapport aux plans initiaux. Ils montrent que cette technique est une technique prometteuse : à ce jour,
nous sommes les premiers à proposer une optimisation continue en les paramètres d’un plan de
traitement en 3D, mais de nombreux développements sont encore nécessaires avant de pouvoir
l’utiliser en clinique.

84
 Chapitre II Optimisation des plans de
traitement en RCMI
Ce chapitre est consacré à l’optimisation des plans de traitement en RCMI qui constitue notre
deuxième axe de travail. L’objectif est de démontrer la faisabilité d'une méthode d'optimisation fondée
sur la technique de calcul de dose rapide que nous avons présentée dans la deuxième partie de ce
mémoire. Cette technique de calcul de dose avec sa reformulation tensorielle nous permet d’aborder le
problème d’optimisation par la méthode de descente de gradient sous contraintes. Dans cette
perspective, nous présenterons tout d’abord le modèle d’optimisation et sa formulation mathématique
dans la section 1. Ensuite, nous allons détailler le calcul du gradient de dose dans la section 2. Enfin,
nous allons montrer dans la section 3 la méthode de descente de gradient utilisée pour résoudre le
problème d’optimisation des plans de traitement.

1 Modèle optimisation et formulation mathématique

Rappelons-nous qu’un plan de traitement est composé de plusieurs faisceaux RCMI positionnés de
façons différentes et que chaque faisceau dispose d’une modulation d’intensité associée. On cherche à
trouver un plan de traitement optimal qui remplit les objectifs cliniques préalablement définis par le
médecin. Idéalement, le système de planification devrait réaliser l’optimisation sur l’ensemble des
paramètres de traitement : la balistique et la modulation. Néanmoins à cause du temps de calcul limité
en clinique, seule la modulation d’intensité de chacun des faisceaux est automatisée. À ce jour,
l’optimisation sur la balistique du traitement est toujours utilisateur-dépendant et doit généralement
être faite manuellement. Il y a parfois une automatisation possible de cette optimisation mais alors
sous la forme d’une sélection de faisceaux parmi un panel très contraint.

Dans le cadre de la thèse, nous proposons une méthode d’optimisation qui optimise non seulement la
modulation d’intensité, mais aussi la balistique du traitement. L’optimisation sur les deux propriétés se
fait en séquentiel ou simultanément. Cette méthode est itérative et déterministe. Elle part d’un plan
issu d’un protocole ou résultant de l’expérience du physicien, qui est considéré comme une bonne
solution a priori. L’objectif sera d’améliorer cette solution potentielle.

Le premier avantage majeur de la méthode est d’être continue en les paramètres avec par défaut tous
les degrés de liberté possible. Ainsi, contrairement aux méthodes proposées dans la littérature, on
explore par défaut un espace de solutions qui permet d’envisager tous les faisceaux possibles pour la
balistique. En d’autres termes, on est capable de placer des faisceaux non-coplanaires et non iso-
centriques. Le second avantage est que l’on peut envisager de limiter la variation de certains
paramètres de façon à respecter des contraintes de valeurs (certains angles impossibles pour l’appareil
clinique), voire à forcer un paramètres à rester constant. Ceci nous permet de spécifier les paramètres à
optimiser en toute liberté. Cette double capacité permet à la méthode d’être exploitable pour n’importe
quel appareillage, y compris comme nous le verrons plus tard pour des accélérateurs hélicoïdaux ou le
Cyberknife.

Nous allons voir dans la section 1.1 la formulation de la fonction d’objectif. Ensuite, les paramètres
d’optimisation pour réaliser cette méthode seront présentés dans la section 1.2 et 1.3.

1.1 Objectifs d’optimisation

Dans le cadre de cette thèse, les objectifs d’optimisation sont fondés sur des critères physiques de
dose: irradier le volume cible conformément à la prescription tout en maintenant la dose aux OARs en
dessous du niveau de tolérance et aussi faible que possible. Ce choix est fait dans le cadre de la

85
faisabilité de dérivation, sachant que les critères biologiques pourraient être aussi utilisés s’ils sont
dérivables (comme l’EUD par exemple).

Les critères de dose : ils comprennent une prescription de dose pour le volume cible et une dose nulle
pour les OAR. Ces prescriptions sont accompagnées par les contraintes de dose qui s’expriment sous
la forme de doses minimales ( ) et maximales ( ) acceptables pour le volume cible et de doses
maximales tolérées pour les OARs.

Ces multiples critères nous permettent de définir un problème d'optimisation multi-objectifs. Pour
optimiser cette fonction, une des démarches courantes consiste à définir une combinaison linéaire de
tous les critères. Pour cela, on doit attribuer un facteur de priorité à chacun d’entre eux. Ces facteurs de
priorité sont divisés par le nombre de points contenus dans chaque structure. Ceci a pour but d’éviter
l’effet dominant de la structure de plus grand volume dans la fonction objectif.

Sous sa forme quadratique, basée uniquement sur la prescription et des contraintes de dose, la fonction
objectif utilisée peut s’écrire :

Où :
est la dose prescrite à la cible,
est la dose de tolérance pour l’OAR,
représente la dose déposée dans le ème voxel du volume cible,
représente la dose déposée dans le ème voxel de l’OAR,
et représentent respectivement les facteurs de priorité attribués à la contrainte cible et à la
contrainte OAR,
et représentent respectivement les nombre de voxels du volume cible et de l’OAR.

L’objectif est de minimiser cette fonction . En effet, pour un plan de traitement donné, la valeur de
quantifie la proximité de la distribution de dose calculée avec la distribution de dose désirée définie
par un ensemble d’objectifs dosimétriques donnés. Il s’agit par conséquent, d’un indicateur chiffré
représentant la qualité d’un plan de traitement, à partir de critères dosimétriques.

1.2 La modulation d’intensité du faisceau

Dans le cadre de cette thèse, on utilise un faisceau RCMI composé de faisceaux élémentaires
de section carrée mis côte à côte selon deux directions. Chaque faisceau élémentaire possède une
fluence uniforme (ce qui revient aussi à parler d’intensité uniforme), illustré sur la Figure 58.

L’ensemble de ces fluences uniformes constitue un vecteur de fluence, noté , pour le faisceau
RCMI. Ce vecteur de fluence doit être positif pour respecter la réalité physique.

86
 Figure 58 La carte de fluence d’un faisceau RCMI

Rappelons-nous que, pour un faisceau RCMI donné, le faisceau élémentaire est indexé par
Donc, la fluence du faisceau élémentaire est désignée par ou 8.

Définir la modulation d’intensité, c’est trouver, pour une balistique définie, le vecteur de fluence qui
minimise la fonction présentée dans la section 1.1.

1.3 Balistique du traitement

Dans le cadre de la thèse, on considère qu’une balistique de traitement est composée de plusieurs
faisceaux RCMI qui sont positionnés dans l’espace trois dimensions. Les grandeurs caractéristiques du
faisceau RCMI sont ici indicées par . Chaque faisceau RCMI est caractérisé dans l’espace par la
position de son origine et par son orientation représentée par ( ; ; ).

Compte tenu que les interactions du faisceau avec l’air sont négligées (assimilé à du vide), le
déplacement du point d’origine selon l’axe du faisceau ne modifie pas le dépôt de dose final9.
Nous limitons donc les déplacements de sur une sphère de rayon qui est suffisamment large
pour envelopper entièrement le fantôme hétérogène. Ceci nous permet de réduire le paramétrage du
déplacement de sur la sphère à deux degré de liberté : (longitude) et (colatitude). et
sont indépendants. est compris entre 0 et et est compris entre 0 et . Le déplacement du
faisceau élémentaire dépend de son , il peut être déduit facilement par une simple
translation de selon l’axe et du faisceau.

L’orientation du faisceau est caractérisée par les angles d'Euler de la variante Tait-Bryan ( ):
, et , qui définissent totalement le triplet ( ; ; ) du faisceau. Les trois angles sont
indépendants. Dans un premier temps, on considère qu’ils sont compris entre – et (Deuxième
partie, Chapitre IV, section 1.3).

Définir une balistique, c’est trouver les paramètres de l’origine ( ; ) et de l’orientation ( ; ;

) de chaque faisceau RCMI indicées par , dans le but de minimiser la fonction présentée dans la
section 1.1. Ces cinq paramètres sont regroupés dans le vecteur paramètre .

2 Gradient de dose
2.1 Composition du gradient de dose
Le gradient est un opérateur qui quantifie au premier ordre la variation d'une fonction induite par les
variations des variables dont elle dépend. Dans notre cas, nous cherchons à connaitre la variation de la
dose déposée en fonction des paramètres du plan de traitement, c’est-à-dire les fluences et les
balistiques.

Étant donnée un plan de traitement composé de faisceaux RCMI, par le principe de superposition,
la dose déposée au point par ces faisceaux RCMI peut s’écrire10 :

8
En pratique, l’ensemble de fluences d’un faisceau RCMI peut s’exprimer sous forme de vecteur ou de matrice
d’où deux façons d’indiquer la fluence d’un faisceau élémentaire.
9
Sur de gros faisceaux, cette assertion est légèrement faux à cause du phénomène de pénombre … mais dans
notre modèle de petits faisceaux mitoyens, ce phénomène de pénombre n’est pas réellement considéré et pris en
compte. De plus, les nouvelles technologies sont d’ailleurs de moins en moins concernées (ex du Cyberknife).
10
La notation est équivalente de la notation pour un point se situe au ième voxel.

87
 ème
, est la dose déposée par le faisceau RCMI

Le gradient de la dose au point par rapport aux paramètres, désigne l’ensemble des paramètres
( ; … … ), des faisceaux RCMI (les paramètres du plan de traitement à optimiser).
Notons que , présente les paramètres du ème faisceau RCMI. est constitué de cinq
paramètres géométriques ( ; ; ; ; ), est le vecteur de fluence composée de
réel positif.

Nous avons :

On sait qu’un faisceau RCMI est composé de faisceaux élémentaires, donc :

où est l’intensité du faisceau élémentaire du ème faisceau RCMI, est la

dose déposée au point par ce faisceau élémentaire avec une fluence unité.

On cherche le gradient de la dose déposée par ce faisceau RCMI. Selon les paramètres de , le
gradient peut s’écrire de deux façons différentes:

Ainsi que,

Ici, on constate que le gradient de la dose par rapport à la matrice de fluence est très facile à calculer, il
suffit de connaitre la dose déposée par chaque faisceau élémentaire. Cette dose peut être calculée par
la méthode Pencil Beam ou Monte Carlo ou autres. C’est pour cette raison que pour une balistique
prédéfinie, la majorité des méthodes d’optimisation utilise la méthode de gradient pour trouver la
matrice de fluence optimale.

Néanmoins, à notre connaissance, il n’existait pas jusqu’à présent de méthode de calcul de dose qui
permette d’expliciter le gradient de la dose par rapport aux paramètres géométriques du faisceau. Par
conséquent, les balistiques du traitement sont définies manuellement, par sélection dans un panel
limité ou par des méthodes purement géométriques mentionnées dans la section 4 du Chapitre I.

En revanche, notre méthode de calcul de dose rapide nous permet d’expliciter un tel gradient par
rapport à la balistique du plan de traitement. Nous allons maintenant continuer la décomposition du
gradient de la dose en tenant compte de la particularité de notre méthode de calcul de dose.

La dose déposée par un faisceau élémentaire, , est calculée par le logiciel Doséclair.
Elle est la somme pondérée de la dose déposée par chaque SBIM :

88
Où est la dose déposée par le SBIM du faisceau élémentaire ( ), est la
pondération du SBIM sur le point .

Dans un premier temps, on considère que la pondération est invariante par rapport aux
changements des paramètres du faisceau. Or ne dépend pas de l’intensité du faisceau.
Nous allons donc nous intéresser uniquement au gradient de la dose déposée par ce faisceau
élémentaire ( ) par rapport aux 5 paramètres géométriques du faisceau RCMI.

La dose déposée par le SBIM est calculée en deux étapes :

D’abord, on calcule la projection du point dans son modèle homogène :

Ensuite, on calcule la dose de cette projection en utilisant un modèle de régression :

Donc, nous avons

Grâce au théorème de dérivation des fonctions composées, le gradient de la dose déposée par le
SBIM peut s’écrire :

Maintenant, nous sommes arrivés à la racine du problème. Pour connaitre le gradient de la dose par
rapport à la balistique du plan de traitement, il suffit de connaitre le gradient du modèle de régression
par rapport au point projeté et le gradient du tenseur de projection par rapport à la balistique du plan de
traitement.

La fonction pour obtenir la dose en milieu homogène finalement utilisée dans ce travail est à base de
splines d’ordre 3 bien adaptées à la nature particulière des courbes à représenter. Elle est donc
modélisée par une grille de polynômes en fonction des coordonnées du point projeté. Pour une

raison de simplification, est exprimé par . Ces polynômes sont de la forme :

. Cette fonction est dérivable (y compris sur les limites de la grille,

par construction des splines) et nous pouvons expliciter ce gradient :

Toutefois si d’autres méthodes de régression que les splines sont employées, il n’y a aucun problème,
la dérivation étant souvent peu courante mais parfaitement possible pour la grande majorité de ces
méthodes (comme les réseaux de neurones artificiels qui avaient aussi été explorés dans les travaux
préliminaires).

89
À ce stade, pour évaluer le gradient de la dose, il ne reste plus qu’un élément à expliciter : le gradient
du tenseur de projection par rapport à la balistique du plan de traitement Nous allons détailler ces
calculs dans les prochaines sections.

2.2 Le gradient du tenseur de projection

La propagation du faisceau RCMI dans le fantôme est réalisée par propagations du faisceau
élémentaire. Chaque propagation du faisceau élémentaire créé ses propres SBIM et donc ses tenseurs
de projections axiales et latérales.

Nous avons vu dans la deuxième partie de ce mémoire que, pour un faisceau élémentaire paramétré
par =( ; ; ; ; ), la projection est représentée par un tenseur d’ordre 2, équivalent à une
matrice, de dimension 4×4. Nous disposons par ailleurs des formules analytique de calculs de ses
tenseurs à partir des cinq paramètres géométriques du faisceau RCMI via des fonctions
trigonométriques. On peut évaluer le gradient de chaque composante de ces tenseurs par rapport à
chaque paramètre en question. Par conséquent, nous pouvons obtenir le gradient de tenseur pour
chaque SBIM et chaque maille voisine.

Les paramètres étant un vecteur ( ; ; ; ; ), et le tenseur de projection étant d’ordre 2, le

gradient du tenseur de projection par rapport aux paramètres est un tenseur d’ordre 3. Toutefois,
pour simplifier les écritures, on va plutôt écrire cinq gradients sous la forme de tenseurs d’ordre 2,
correspondant au gradient par rapport à un seul des paramètres.

Ainsi, le gradient du tenseur par rapport à , noté , est appelé un tenseur de gradient . Ce dernier
nous permet de connaitre la pente de la variation des coordonnées du point projeté dans le modèle
homogène si on fait varier l’angle du faisceau RCMI. Évidemment, cette proposition est aussi
valable pour les autres paramètres. Donc, pour un SBIM (ou une maille voisine), nous cherchons à
exprimer les cinq tenseurs de gradient associés.

Dans cette perspective, nous avons commencé par exprimer les cinq tenseurs de gradient de la maille
d’entrée. Ensuite en suivant le même principe de propagation, nous avons établis le mécanisme de
propagation en axial et en latéral. Le calcul est détaillé en Annexe C et Annexe D.

Finalement, on obtient un ensemble de tenseurs de gradient analytiques en fonction des paramètres

géométriques du faisceau RCMI. En tenant compte de la composition du gradient de dose présenté
dans la section précédent, nous pouvons connaitre l’importance de la variation de la dose en chaque
point du fantôme en faisant varier la balistique du plan de traitement. Cette information précieuse nous
donne la possibilité de résoudre le problème d’optimisation de la balistique par la méthode de descente
de gradient.

3 Descente de gradient
3.1 Algorithme générale
La technique de recherche de la solution optimale par la méthode du gradient est la plus couramment
utilisée en radiothérapie car une de ces principales particularités est sa rapidité d’obtention d’une
solution acceptable.

90
Le processus d’optimisation va chercher à minimiser la valeur de la fonction objectif , ce qui est
souvent fait en recherchant les valeurs de qui annule la dérivée11 (gradient au cas des variables
multiples).
min

Dans le cas de la thèse, on définit une fonction objectif quadratique qui mesure l’écart entre la dose
déposée et la dose prescrite. Dans le cas le plus simple, la dose prescrite pour les OARs est nul. La
variable est constituée par l’ensemble des paramètres géométrique des faisceaux et leurs intensités
associées.

La recherche de se fera de façon itérative, selon la relation générale des gradients :

= +

où désigne le numéro de l’itération, présente la taille du pas de recherche (variable ou non) et

est la direction de descente.

Comme notre objectif n’est pas de créer un nouvel algorithme d’optimisation, nous avons fait le choix
d’utiliser pour notre problématique, la méthode des gradients dite de plus forte pente (« steepest
descent »). Comme son nom indique, cette méthode utilise la direction de descente de plus forte pente
afin d’atteindre un minimum le plus rapidement possible, soit :

D’où :

Nous avons utilisé un pas fixe dans les simulations présentées dans ce mémoire.

Comme présenté précédemment dans la section 2, à chaque point du fantôme, nous disposons du
gradient de dose par rapport aux paramètres du plan de traitement. Nous pouvons optimiser la fonction
par rapport à chaque paramètre avec la contrainte des intensités des faisceaux non négatives. Cette
contrainte est respectée par une remise à zéro forcée lorsque l’intensité se trouve négative pendant le
processus.

Notons que notre but est de trouver une solution optimale à partir d’un point de départ pas trop mal
initialisé. Donc, le choix du n’est pas arbitraire. Ensuite, au fur et à mesure des itérations, la valeur
de la fonction objectif décroît pour atteindre le minimum. Nous arrêtons le processus lorsque le
gradient est suffisamment petit ou que le nombre maximum d’itérations est atteint.

11
Ce qui si l’on considère l’ensemble des fonctions mathématiques n’est pas strictement équivalent puisque ce
n’est qu’un critère d’extrémum : le gradient est également nul pour les maximums, les plateaux et les points-
selle. Toutefois la nature particulière des fonction objectifs (souvent des sommes à coefficients strictement
positifs de fonctions quadratiques comme pour ) et le principe même de la méthode (descente de gradient) font
qu’effectivement, sous ces hypothèses, atteindre un point de gradient nul est équivalent à atteindre un minimum
local.

91
 Chapitre III Implémentation logicielle et
validation des résultats
Dans ce chapitre, nous présenterons tout d’abord l’intégration des calculs de gradient dans le prototype
existant. Nous verrons ensuite les évaluations des paramètres individuelles pour valider les gradients
implémentés. Nous détaillerons enfin les résultats d’optimisation concernant d’une part, un fantôme
homogène d’eau avec cible en forme de C et d’autre part, un fantôme hétérogène poumon-os avec une
tumeur plus simple. Pour le fantôme homogène, la comparaison des résultats obtenus par des faisceaux
RCMI et de type Cyberknife, sera présentée, montrant la flexibilité des méthodes.

1 Implémentation logicielle
L’implémentation de la méthode présentée dans cette thèse s’est faite en deux parties : simulation et
optimisation. D’abord, pour pouvoir être efficace au niveau de la rapidité de simulation, la méthode de
calcul de dose a été programmée en langage C++. Ce prototype s’appelle ‘Doséclair’. Il s’agit d’un
exécutable s’exécutant en ligne de commande. Les résultats des calculs sont sauvegardés dans un
fichier spécifique. Les tâches effectuées par ce prototype sont les suivantes:

- Charger un fantôme et définir sa composition de matériaux

- Charger les modèles homogènes
- Définir les faisceaux
- Faire la simulation
- Sauvegarder les résultats

J’ai travaillé à partir du code de Baptiste Blanpain, en intégrant les nouvelles fonctionnalités
développées dans le cadre de cette thèse. Pour ce faire, j’ai implémenté une nouvelle classe ‘Tenseur’
qui a pour but d’automatiser les calculs de fonction de projection par les calculs tensoriels et en même
temps calculer les gradients associés. Cette implémentation est intégrée dans les codes qui gèrent la
propagation du faisceau dans le fantôme. À cet effet, la propagation a nécessité également des
modifications et des adaptations importantes.

En outre, comme nous l’avons vu dans le précédent chapitre, un atout majeur de notre méthode est de
pouvoir calculer la dose en certains points du fantôme seulement. De ce fait, j’ai réalisé un module
pour pouvoir retourner non pas la dose dans le fantôme entier mais la dose en des points spécifiés du
fantôme. Cette fonctionnalité est utilisée dans la phase d’optimisation pour accélérer les calculs.

Une fois les routines de simulation complètes, nous avons cherché à créer un module d’optimisation
destiné au plan de traitement. Pour avoir une grande souplesse d’analyse a posteriori et une bonne
visualisation des résultats, nous avons choisi d’implémenter ce module sous MATLAB.12

Un fichier MEX est utilisé comme une interface entre les routines de simulation C++ et le module
d’optimisation MATLAB. Comme une fonction MATLAB ordinaire, on l’appelle en lui passant des
arguments d’entrée, qui sont, dans notre cas :

- Les paramètres géométriques du fantôme

12
Pour pouvoir utiliser les routines C++ existantes, j’ai couplé le code C++ de Doséclair avec
MATLAB par l’intermédiaire d’un fichier MEX, qui est l’acronyme de MATLAB Executable. Ce
fichier MEX est écrit en C++, qui une fois compilé, est utilisé dans MATLAB de la même manière
que les fichiers M.

92
 - La composition des matériaux du fantôme
- Les paramètres du faisceau RCMI

Et en récupérant des arguments de sortie, qui sont, dans notre cas :

- La distribution de dose calculée par Doséclair

- Les gradients de dose associés

Dans le module d’optimisation MATLAB, cette fonction MEX est appelée à chaque itération de la
boucle d’optimisation, elle prend en entrée les paramètres à optimiser (paramètres du faisceau RCMI)
et retourne la distribution de dose et son gradient associé. Ces derniers vont être utilisés pour calculer
la fonction objectif et la mise à jour des paramètres.

2 Validation des paramètres individuels

Dans le but de valider les calculs de gradients présentés dans le Chapitre II, nous avons étudié
individuellement l’évolution des deux groupes de paramètres géométriques représentatifs du plan de
traitement pendant l’optimisation. Le premier groupe contient les paramètres et , qui définissent le
positionnent du point d’origine du faisceau, le deuxième groupe contient les paramètres , et , qui
définissent l’incidence du faisceau. La matrice de fluence sera validée dans la section 3.3.1
séparément.

La variation de chaque paramètre géométrique apporte une modification de la balistique

indépendamment de l’un à l’autre. Cette modification de la balistique entraine une modification de la
distribution de dose que, théoriquement, nous sommes capables de prédire. Pour vérifier que cette
prédiction est correctement modélisée à l’aide de la reformulation tensorielle, une série de test est
réalisé. Chaque test est désigné pour valider un seul paramètre sauf dans le cas de l’origine du faisceau
défini par le couple ( , ). Chaque test est initialisé de telle sorte que le faisceau ne passe pas dans le
volume cible. Le but est de savoir si en modifiant un seul paramètre (ou ( , ) simultanément) selon
une descente de gradient simple, la balistique du faisceau se stabilise par déposer une forte dose dans
le volume cible. La localisation du volume cible est différente en fonction du paramètre à valider.
Nous considérons que, pour un paramètre donné, un test validé implique une bonne modélisation du
gradient par rapport à ce paramètre.

2.1 Conditions de simulation

Le fantôme utilisé dans cette section est composé d’une seule maille homogène d’eau de 12×12×12
cm3 (60×60×60 voxels). L’origine du repère fantôme se situe au centre de ce fantôme.

Un faisceau RCMI (noté ) composé de 4 faisceaux (de section carré de coté 1 cm) mis en carré
est utilisé.

Pour la validation des paramètres géométriques, le faisceau élémentaire , marqué par un carré
rose sur la Figure 59, est utilisé. En pratique, ce faisceau élémentaire a une fluence unitaire tandis que
les trois autres faisceaux ont une fluence nulle.

Figure 59 Faisceau RCMI et faisceau élémentaire utilisés pour la validation

93
Dans cette section, le volume cible est simplifié en une ‘zone cible’ : une zone 2D située sur le plan de
coupe qui est à 2cm de profondeur dans le fantôme. La fonction de coût est quadratique sous la
forme :

Où :
est la dose prescrite à la cible,
représente la dose déposée dans le ième voxel de la zone cible.

2.2 Validation du positionnement du point d’origine du faisceau

Nous avons vu dans le précédent chapitre que l’origine du faisceau RCMI est positionnée sur une
sphère qui enveloppe le fantôme. L’objectif de cette simulation est d’observer le déplacement du point
sachant que l’on veut déposer la dose prescrite dans la zone cible, les autres paramètres restant
invariants pendant l’optimisation.

À l’initialisation, le faisceau RCMI est placé de telle sorte que son faisceau élémentaire rentre
orthogonalement à la surface du fantôme et passe par son centre. La courbe isodose13 initiale est
centrée sur le plan de coupe et est en dehors de la zone cible encadrée par un rectangle pointillé
noir, illustré en Figure 60(a). Sur la Figure 60(b), (c) et (d), nous pouvons constater l’évolution de
cette courbe isodose, au fur à mesure des itérations, elle se rapproche de plus en plus à la zone cible et
enfin rentre dans cette dernière. Ceci est quantifié par une fonction objectif décroissante qui converge
vers 0.0025, illustrée en Figure 61(a). Sur la Figure 61(b) et (c), nous pouvons remarquer que le
paramètre partant de 175.2° est descendu à 166.5°, le paramètre , au contraire, n’a pas eu un
changement significatif, à la convergence, reste à 45° à trois décimal près. Si on regarde l’évolution
de l’abscisse et de l’ordonnée du point d’origine du faisceau , ils sont initialisé à (0, 0) et se
retrouve en (0.9055, 0.9056). Ceci est cohérent avec un déplacement du volume traité dans la zone
cible, qui se situe à environ 1 cm plus haut et à droite du centre du fantôme.

Les vues 3D avant et après l’optimisation sont présentées en Figure 62 et Figure 63 respectivement. Le
carrée rouge représente la zone cible, le bâton en bleu foncé représente le faisceau , la maille
d’eau est représenté par une cube en bleu transparent.
6 6

3 3
X (cm)

X (cm)

0 0

-3 -3

-3 0 3 6 -3 0 3 6
Y (cm) Y (cm)

(a) (b)

13
Dans cette section, la courbe isodose est définie par isodose 100%, c'est-à-dire qu’elle entoure la zone où la
dose prescrite est au moins atteinte.

94
 6 6

X (cm) 3 3

X (cm)
0 0

-3 -3

-3 0 3 6 -3 0 3 6
Y (cm) Y (cm)

(c) (d)

Figure 60 Validation du point d’origine : courbe isodose initiale (a) et à l’itération 10 (b), 30 (c), et 60 (d)

Fonction objectif Evolution du paramètre Phi Evolution du paramètre Theta

0.09 176 45.0025

0.08 175

174 45.002
0.07

173
0.06

172 45.0015
0.05

171

0.04
170 45.001

0.03
169

0.02
168 45.0005

0.01 167

0 166 45
0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60
Nombre d'itération Nombre d'itération Nombre d'itération

(a) (b) (c)

Figure 61 Fonction objectif (a) et les évolutions des deux paramètres ϕ (b) et θ (c)

(a) (b)

95
 (c) (d)

Figure 62 Vues initiales de la validation du point d’origine: vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

(a) (b)

(c) (d)

Figure 63 Vues finales de la validation du point d’origine : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

2.3 Validation de l’orientation du faisceau

96
Les angles d'Euler sont utilisés pour caractériser l’orientation du faisceau. Ils sont présentés dans la
deuxième partie, Chapitre IV, section 1.3 du mémoire. Nous avons choisi une variante (nommée Tait-
Bryan z-y-x). Dans la variante choisie, la succession de rotations en défini par :

- une rotation par autour de l'axe ,

- une rotation par autour du nouvel axe (parfois aussi noté )
- une rotation par autour du nouvel axe (parfois aussi noté ).

2.3.1 Paramètre

Le paramètre désigne l’angle de rotation du repère faisceau autour de l’axe du repère fantôme. En
partant du même principe, l’objectif de cette simulation est d’observer l’évolution du paramètre pour
déposer la dose dans la zone cible. Dans ce cas, on empêche les variations des autres paramètres
pendant l’optimisation.

À l’initialisation, le faisceau RCMI est placé de la même manière que la précédente simulation. C’est-
à-dire le faisceau élémentaire rentre orthogonalement à la surface du fantôme et passe par son
centre. En revanche, la zone cible est placé à 1 cm plus haut selon l’axe par rapport au centre du plan
de coupe , illustré en Figure 64(a)14. Plus précisément, le centre de la zone cible se situe au (1, 0) sur
ce plan de coupe. Sur la Figure 64(b), (c) et (d), nous pouvons constater que la courbe isodose se
rapproche à la zone cible pendant les itérations et à la fin elle rentre dans la zone. Nous pouvons
constater aussi que la fonction objectif est décroissante et qu’elle converge au bout de 30 itérations,
illustré en Figure 65(a). Sur la Figure 65(b), le paramètre converge vers -90°. Dans la réalité, = -
90° représente une rotation de 90° du repère faisceau ( ; ; ) autour de l’axe dans le sens
inverse des aiguilles d’une montre (voir Figure 66). Cette rotation est équivalente à un déplacement
d’1 cm du point d’origine du faisceau selon l’axe , ce qui déplace là où était et explique
le dépôt de dose dans la zone cible (cf. Figure 66).

Les vues 3D avant et après l’optimisation sont présentées en Figure 67 et Figure 68 respectivement.

6 6

3 3
X (cm)

X (cm)

0 0

-3 -3

-3 0 3 6 -3 0 3 6
Y (cm) Y (cm)

(a) (b)

14
La zone définie à la simulation précédente n’est en effet pas atteignable par une simple rotation de , alors que
cette nouvelle zone oui.

97
 6 6

3 3
X (cm)

X (cm)
0 0

-3 -3

-3 0 3 6 -3 0 3 6
Y (cm) Y (cm)

(c) (d)

Figure 64 Validation du paramètre : courbe isodose initiale (a) et à l’itération 10 (b), 30 (c), et 50 (d)

Fonction objectif Evolution du paramètre Alpha

0.05 0

0.045 -10

0.04
-20

0.035
-30

0.03
-40

0.025

-50
0.02

-60
0.015

-70
0.01

-80
0.005

0 -90
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Nombre d'itération Nombre d'itération

(a) (b)

Figure 65 Fonction objectif (a) et l’évolution du paramètre (b)

Figure 66 Le déplacement du point d’origine avec l’évolution du paramètre

98
 (a) (b)

(c) (d)

Figure 67 Vues initiales de la validation du paramètre : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

(a) (b)

99
 (c) (d)

Figure 68 Vues finales de la validation du paramètre : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

2.3.2 Paramètre

Le paramètre désigne l’angle de rotation du repère faisceau autour de l’axe du repère fantôme15.
L’objectif de cette simulation est d’observer l’évolution du paramètre pour déposer la dose dans la
zone cible. Dans ce cas, on empêche les variations des autres paramètres pendant l’optimisation.

Le faisceau RCMI est initialisé de la même manière que la précédente simulation. La zone cible est
placé à environs 1 cm plus bas selon l’axe par rapport au centre du plan de coupe , illustré en
Figure 69(a). Plus précisément, le centre de la zone cible se situe au (-1, 0) sur ce plan de coup. Nous
pouvons constater que la fonction objectif est décroissante et elle converge au bout de 5 itérations,
illustré en Figure 70(a). La Figure 69(b) montre que la courbe isodose est rentrée dans la zone cible
après la convergence de la fonction objectif. Cette convergence est réalisée grâce à une évolution du
paramètre qui tend vers 10°. Sur la Figure 71, on constate que = 10° représente une rotation de
10° du repère faisceau ( ; ; ) autour de l’axe dans le sens des aiguilles d’une montre, le trait
en pointillé noir représente la localisation de la zone cible sur le plan de coupe , cette rotation permet
de faire incliner le faisceau de sorte qu’il passe dans la zone cible, ce qui est le résultat attendu
de cette simulation.

Les vues 3D avant et après l’optimisation sont présentées en Figure 72 et Figure 73 respectivement.

15
En théorie, c’est , le nouvel axe après la première rotation autour de . Ici la rotation autour de l’axe est
nulle, l’axe est équivalent à l’axe .

100
 6 6

3 3
X (cm)

X (cm)
0 0

-3 -3

-3 0 3 6 -3 0 3 6
Y(cm) Y (cm)

(a) (b)

Figure 69 Validation du paramètre : courbe isodose initiale (a) et à l’itération 10 (b)

Fonction objectif Evolution du paramètre Beta

0.04 10

9
0.035

8
0.03
7

0.025
6

0.02 5

4
0.015

3
0.01
2

0.005
1

0 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Nombre d'itération Nombre d'itération

(a) (b)

Figure 70 Fonction objectif (a) et l’évolution du paramètre (b)

Figure 71 Configuration géométrique avec paramètre = 0° et = 10°

101
 (a) (b)

(c) (d)

Figure 72 Vues initiales de la validation du paramètre : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

(a) (b)

102
 (c) (d)

Figure 73 Vues finales de la validation du paramètre : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

2.3.3 Paramètre

Le paramètre désigne l’angle de rotation du repère faisceau autour de l’axe du repère fantôme16.
L’objectif de cette simulation est d’observer l’évolution du paramètre pour déposer la dose dans la
zone cible. Dans ce cas, on empêche les variations des autres paramètres pendant l’optimisation.

Même principe d’initialisation que la précédente simulation. La zone cible est placé à 1 cm à droite
selon l’axe par rapport au centre du plan de coupe . Plus précisément, le centre de la zone cible se
situe au (0, 1) sur ce plan de coupe. Nous pouvons constater que la fonction objectif est décroissante et
elle converge au bout de 10 itérations, illustré en Figure 74(a). La Figure 75(b) montre que la courbe
isodose est rentrée dans la zone cible après la convergence de la fonction objectif. Cette fois ci, la
convergence est réalisée grâce à une évolution du paramètre qui tend vers 10°. Sur la Figure 76, on
voit que = 10° représente une rotation de 10° du repère faisceau ( ; ; ) autour de l’axe dans
le sens inverse des aiguilles d’une montre, cette rotation permet de faire incliner le faisceau de
sorte qu’il passe par le trait en pointillé noir représentant la zone cible. L’objectif de cette simulation
est atteint.

Les vues 3D avant et après l’optimisation sont présentées en Figure 77 et Figure 78 respectivement.

6 6

3 3
X (cm)

X (cm)

0 0

-3 -3

-3 0 3 6 -3 0 3 6
Y (cm) Y (cm)

(a) (b)

16
En théorie, c’est , le nouvel axe après double rotation. Ici ils sont identiques car les deux rotations
précédentes sont nulles.

103
 Figure 74 Validation du paramètre : courbe isodose initiale (a) et à l’itération 30(b)

Fonction objectif Evolution du paramètre Gamma

0.04 10

9
0.035

8
0.03
7

0.025
6

0.02 5

4
0.015

3
0.01
2

0.005
1

0 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Nombre d'itération Nombre d'itération

(a) (b)

Figure 75 Fonction objectif (a) et l’évolution du paramètre (b)

Figure 76 Configuration géométrique avec paramètre = 0° et = 10°

(a) (b)

104
 (c) (d)

Figure 77 Vues initiales de la validation du paramètre : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

(a) (b)

(c) (d)

Figure 78 Vues finales de la validation du paramètre : vue de dessus (a), vue de côté (b), vue de face (c), vue
isométrique (d)

3 Fantôme homogène d’eau avec cible de forme concave

3.1 Conditions de simulation

105
Nous avons validé nos gradients sur les simulations simplifiées dans la section précédente, maintenant,
nous allons étudier un cas plus complexe.

Dans cette section, nous allons conserver le fantôme homogène d’eau utilisé dans les précédentes
simulations. Nous choisissons cette fois ci un volume cible 3D de forme concave autour d’un organe à
risque. Cette configuration est représentative de cas couramment traités en clinique même si elle est
simplifiée ici : un cancer de la prostate à proximité du rectum, le lymphome à proximité de la moelle
épinière ou bien le méningiome à proximité des yeux.

Le volume cible – ayant une forme en ‘C’ – enveloppe un organe à risque sous forme d’un
parallélépipède rectangle. L’ensemble se situe au centre du fantôme. La fonction de coût est
quadratique de la forme :

Où :
est la dose prescrite à la cible,
représente la dose déposée dans le ième voxel du volume cible,
représente la dose déposée dans le jème voxel de l’OAR,
et représentent respectivement les nombre de voxels du volume cible et de l’OAR.

Nous allons dérouler nos simulations sur deux types de traitement : le Cyberknife et RCMI classique.

3.2 Cyberknife - faisceaux élémentaires indépendants

Le traitement du type Cyberknife est caractérisé par sa grande flexibilité de ciblage et sa souplesse
dans l’orientation des faisceaux. Il permet de réaliser aisément des traitements avec des multiples
faisceaux non coplanaires et non isocentriques qui convergent vers le volume cible. Les collimateurs
circulaires de diamètre différents sont utilisés pour définir la forme du faisceau.

Or un des avantages de notre méthode est de pouvoir optimiser des faisceaux avec une très grande
souplesse dans la balistique du traitement, nous allons d’abord étudier cette simulation avec des
faisceaux du type Cyberknife (faisceau élémentaire). Il faut noter que, à ce jour, nous n’avons pas
encore implémenté le faisceau circulaire dans Doséclair, nous allons utiliser les faisceaux de section
carrée pour la simulation.

3.2.1 Intensité uniforme (49 mins)

Pour commencer, nous allons étudier l’évolution des paramètres géométrique du faisceau, c'est-à-dire
que chaque faisceau a une fluence uniforme et invariante pendant l’optimisation.

Dans cette simulation, 10 faisceaux élémentaires non coplanaires sont utilisés. Chaque faisceau
possède un repère faisceau défini par son point d’origine (le point focal de l’accélérateur) et les
trois axes du repère ( ; ; ) dont le vecteur désigne la direction d’irradiation du faisceau.
Un total de 50 paramètres est à optimiser simultanément. La dose prescrite est de 50Gy, le nombre
maximum d’itération est fixé à 600. Nous avons utilisé un pas fixe dans cette simulation.

À l’initialisation, 5 faisceaux coplanaires et isocentriques sont placés dans le plan Y – Z et X –Y

respectivement. Ces 5 faisceaux se positionne sur 1, 5, 7, 9, 11 heure d’une horloge, schématisé en
Figure 79. Les deux plans passent par le centre du fantôme. Les rectangles bleus composent la coupe
du volume cible et le carré rose représente la coupe de l’OAR. Les dix axes des faisceaux sont
présentés par les flèches grises et ils sont dans ces deux plans respectivement. Le cadre rose de la
Figure 79 (b) représente la projection de l’OAR sur ce plan qui néanmoins n’intersecte pas l’OAR.

106
 Vue Y - Z Vue X - Y

(a) (b)
Figure 79 La balistique initiale dans le plan de coup Y –Z (a) et X – Y (b) : fantôme homogène d’eau avec la cible de
forme ‘C’, faisceau type Cyberknife sans variation d’intensité

En partant de cette initialisation, la Figure 80 montre l’évolution des distributions de dose (illustré via
des plans de coupe et des courbes d’isodoses) pendant le processus d’optimisation. Dans le cadre de
cette étude, la fonction d’objectif ne prend en compte que le volume cible et l’OAR, nous avons choisi
de ne calculer la dose que dans les deux entités anatomiques.

Sur la Figure 80, le plan de coupe est défini par l’axe et du repère fantôme. Sur chaque plan de
coupe, nous pouvons voir que le volume cible en forme de ‘C’ est composé par trois zones
rectangulaires encadrées en pointillé blanc, l’OAR est encadré en pointillé rouge. La dose n’est
calculée que dans ces deux régions, c’est la raison pour laquelle le reste du plan est coloré en bleu
foncé qui normalement signifie une dose nulle.

Sur la Figure 80(a), nous pouvons constater que, à l’initialisation, les faisceaux déposent une très forte
dose (surdosage) qui s’élève à 100Gy, présenté par le chiffre 0.1 sur l’échelle de couleur, sur deux
zones symétriques dans le volume cible. Les autres zones de la cible ont subi un sous-dosage. Au bout
de la 100ème itération, illustré par Figure 80(b), la dose est mieux répartie dans le volume cible et
l’OAR est mieux épargné, illustré par Figure 80(b bis). Au bout de la 300ème itération (Figure 80(c)), la
dose est déposée avantageusement dans la zone rectangulaire à gauche de la cible sachant que cette
zone était en sous-dosage à l’initialisation. Malgré cela, quelques surdosages présents dans la cible et
l’OAR est un peu plus irradié. Au bout de la 600ème itération, ces surdosages dans la cible ont été
résolus sans avoir dégradé le dosage dans l’OAR.

107
 Plan de coup Y - Z à l'initialisation Isodose à l'initialisation
0.1

0.09
-4 -4
0.08

-2 -2 0.07

0.06

Z (cm)
Z (cm)

0 0 0.05

0.04
2 2
0.03

4 4 0.02

0.01
6 6
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)

(a) (a bis)

Plan de coup Y - Z à 100 ème itération Isodose à 100 ème itération

0.1

0.09
-4 -4
0.08

-2 -2 0.07

0.06
Z (cm)

Z (cm)

0 0
0.05

0.04
2 2
0.03

4 4 0.02

0.01
6 6
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)
(b) (b bis)

Plan de coup Y -Z à 300 ème itération Isodose à 300 ème itération

0.1

0.09
-4 -4
0.08

-2 -2 0.07

0.06
Z (cm)

Z (cm)

0 0 0.05

0.04
2 2
0.03

4 4 0.02

0.01
6 6
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)

(c) (c bis)

108
 Plan de coup Y - Z à 600 ème itération Isodose à 600 ème itération
0.1

0.09
-4 -4
0.08

-2 -2 0.07

0.06
Z (cm)

Z (cm)
0 0 0.05

0.04
2 2
0.03

4 4 0.02

0.01
6 6
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)
(d) (d bis)

Figure 80 Les plans de coupe et les courbes isodoses : fantôme homogène d’eau avec la cible de forme ‘C’, faisceau
type Cyberknife sans variation d’intensité

Maintenant, nous allons comparer les deux histogrammes dose-volume avant et après l’optimisation,
(Figure 81).

Tout d’abord, nous pouvons constater que le pourcentage de volume cible qui reçoit une dose
supérieure à 50Gy (la dose prescrite) est passée de 30% à 40%. Le volume recevant une dose
supérieure à 30Gy (60% de la dose prescrite) est passé de 40% à 70%. De même, le pourcentage de
volume a bien augmenté pour chaque niveau de dose. Ceci montre que l’optimisation a mieux
homogénéisé la dose déposée dans la cible. Certes, seulement 40% du volume de la cible ayant la dose
prescrite n’est pas suffisant pour un plan final. En réalité, un plan bien initialisé avec un nombre de
faisceau suffisant et bien positionné est nécessaire pour avoir un plan final satisfaisant. Le choix d’une
fonction objectif et d’une méthode de descente de gradient plus sophistiquée peuvent aussi influencer
le résultat final. Dans le cadre de la thèse, nous ne cherchons pas à trouver un plan final utilisable en
clinique, mais nous cherchons plutôt à améliorer un plan initial.

Deuxièmement, nous pouvons constater que l’OAR est mieux épargné à la fin de l’optimisation. La
dose maximale reçue est descendue à 23Gy. Le pourcentage de volume d’OAR qui reçoit une dose
supérieure à 10Gy (20% de la dose prescrite) est passé de 48% à 15%.

Ces améliorations sont confirmées par l’évolution de la fonction objectif (Figure 82). La fonction est
décroissante et elle a subi une forte descente pendant les 300 premières itérations et converge
lentement par la suite.

En conclusion, le plan de traitement peut être amélioré par rapport au plan initial lorsqu’on fait varier
les incidences des faisceaux à petit degré en continu. Nous allons voir l’influence de l’intensité du
faisceau sur cette simulation dans la prochaine section.

109
 Histogramme dose-volume Histogramme dose-volume
1 1
Cible Cible
0.9 OAR 0.9 OAR

0.8 0.8

0.7 0.7
Volume [%]

Volume [%]
0.6 0.6

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.01 0.02 0.03 0.04 0.05 0.06
Dose Dose
(a) (b)

Figure 81 Les histogrammes dose-volume initiale (a) et finale (b) : fantôme homogène d’eau avec la cible de forme ‘C’,
faisceau type Cyberknife sans variation d’intensité

x 10
-3 Fonction objectif
9

1
0 100 200 300 400 500 600
Nombre d'iteration

Figure 82 La fonction objectif : fantôme homogène d’eau avec la cible de forme ‘C’, faisceau type Cyberknife sans
variation d’intensité

3.2.2 Intensité variables (24.5 mins)

Dans cette section, nous allons optimiser le même plan mais avec une initialisation différente, illustré
en Figure 83. Les faisceaux traversent quasi-orthogonale le fantôme. Chaque carré gris représente la
projection géométrique du faisceau dans le plan de coup. Le gris claire signifie que le faisceau propage
vers le lecteur. Le gris foncé signifie que le faisceau propage en s’éloignant le lecteur. Un biais est
introduit pour chaque incidence pour que les faisceaux ne soient pas parfaitement orthogonaux à
l’entrée du fantôme à l’initialisation.

Notons que la fluence de chaque faisceau est aussi à optimiser pendant le processus, ce qui fait une
optimisation simultanée de 60 paramètres. Le nombre maximum d’itérations est fixé à 300. Nous
avons utilisé un pas fixe pour les paramètres géométrique et la fluence dans cette simulation.

110
 Vue Y - Z Vue X - Y

(b)

(b)

(a) (b)
Figure 83 La balistique initiale dans le plan de coup Y –Z (a) et X – Y (b) : fantôme homogène d’eau avec la cible de
forme ‘C’, faisceau type Cyberknife avec variation d’intensité

On affiche le plan de coupe défini par l’axe et du repère fantôme et ses courbes d’isodoses
associées. Sur chaque plan de coupe, le volume cible est marqué par deux zones rectangulaires
encadrées en pointillé blanc, l’OAR est encadré en pointillé rouge.

Plan de coup X - Y à l'initialisation Isodose à l'initialisation

0.06

-4 -4
0.05

-2 -2 0.04
X (cm)

X (cm)

0 0 0.03

2 2 0.02

4 4 0.01

-4 6 0
-4 -2 0 2 4 -4 -4 -2 0 2 4 6
Y (cm) Y (cm)

(a) (a bis)

Plan de coup X - Y à 30ème itération Isodose à 30ème itération

0.06

-4 -4
0.05

-2 -2 0.04
X (cm)

X (cm)

0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)

(b) (b bis)

111
 Plan de coup X - Y à 100ème itération Isodose à 100ème itération
0.06

-4 -4
0.05

-2 -2 0.04
X (cm)

X (cm)
0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)

(c) (c bis)

Plan de coup X - Y à 300ème itération Isodose à 300ème itération

0.06

-4 -4
0.05

-2 -2 0.04
X (cm)
X (cm)

0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)

(d) (d bis)

Figure 84 Les plans de coup et les courbes isodose : fantôme homogène d’eau avec la cible de forme ‘C’, faisceau type
Cyberknife avec variation d’intensité

Rappelons-nous que la dose prescrite pour la cible est 50 Gy (présenté par le chiffre 0.05 sur l’échelle
de couleur). Sur la Figure 84(a), nous pouvons constater que, à l’initialisation, les faisceaux déposent
une très faible dose dans le volume cible, en particulier, la dose dans la partie haute de la cible à droite
est quasi nulle. Au bout de la 30ème itération (Figure 84 (b)), la dose augmente dans le volume cible et
l’OAR est mieux épargné (Figure 84 (b bis)). Au bout de la 100ème itération (Figure 84 (c)), la dose
continue à augmenter dans la cible, et elle commence à être homogénéisée dans la cible notamment à
droite dans laquelle il existait un sous dosage au début. Ce « remplissage » de dose est réalisé grâce
aux évolutions des paramètres géométriques qui réorientent les faisceaux. Au bout de la 300 ème
itération (Figure 84 (d)), la dose atteint 50 Gy dans la majorité du volume cible et l’OAR est mieux
épargné par rapport à l’initialisation. L’augmentation de la dose dans le volume cible est obtenue,
principalement, à l’aide de l’évolution de fluence de chaque faisceau. Cette dernière est passée de 0.2
pour tous les faisceaux à l’initialisation à une valeur entre 0.7 et 0.9 selon les faisceaux en fin
d’optimisation.

112
 Histogramme dose-volume Histogramme dose-volume
1 1
Cible Cible
0.9 OAR 0.9 OAR

0.8 0.8

0.7 0.7
Volume [%]

Volume [%]
0.6 0.6

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.01 0.02 0.03 0.04 0.05 0.06
Dose Dose
(a) (b)

Figure 85 Les histogrammes dose-volume initiale (a) et finale (b) : fantôme homogène d’eau avec la cible de forme ‘C’,
faisceau type Cyberknife avec variation d’intensité

Si on regarde les deux histogrammes dose-volume, illustrés en Figure 85, nous pouvons constater que
la totalité du volume cible a reçu une dose supérieure à 20 Gy sachant que cette dose n’est présente
que dans 10% du volume cible à l’initialisation. En outre, 45% du volume cible a reçu la dose prescrite
(50 Gy) et plus de 80% de volume a reçu une dose supérieur à 40 Gy en fin d’optimisation. Ces
niveaux de dose ne sont présents nul part dans le volume cible à l’initialisation.

Pour l’OAR, la dose maximale est réduite à 5 Gy après l’optimisation sachant qu’il y a 16% du
volume d’OAR qui recevait une dose supérieure à 10Gy à l’initialisation.

Si on compare le plan obtenu dans la section précédente avec ce dernier, on peut conclure que ce
dernier entraine une augmentation significative de l’homogénéité au sein du volume cible et une forte
diminution en absolu de la dose reçue par l’OAR.

Ces observations sont confirmées par la fonction objectif ci-dessous (Figure 86). La fonction converge
vers une valeur quasi nulle au bout de 300 itérations tandis que cette valeur est 6 fois plus grande dans
la section précédente. Entre autres, alors qu’il y a beaucoup plus de paramètres à optimiser, le calcul
est accéléré grâce à la diminution du nombre d’itérations nécessaires.

x 10
-3 Fonction objectif
4.5

3.5

2.5

1.5

0.5

0
0 50 100 150 200 250 300
Nombre d'iteration
Figure 86 La fonction objectif : fantôme homogène d’eau avec la cible de forme ‘C’, faisceau type Cyberknife avec
variation d’intensité

En conclusion, une bonne initialisation du plan joue un rôle important pour l’optimisation. La
participation de la fluence dans l’optimisation augmente les degrés de liberté du problème et elle
apporte un bénéfice pour le plan final au niveau de la dose dans la cible.

113
3.3 Faisceaux RCMI classique
Dans cette section, nous avons utilisé 6 faisceaux RCMI de section 4cm×4cm (16 faisceaux
élémentaires par chaque faisceau RCMI). L’objectif est d’observer d’une part l’évolution des
incidences des 6 faisceaux, d’autre part l’évolution de la matrice de fluence de chaque faisceau RCMI.

3.3.1 Plan angles fixes avec modulation (3.7 mins)

Nous étudions d’abord le cas le plus simple qui est mis en œuvre en fixant les incidences des faisceaux
RCMI. Ces incidences sont définies de telle sorte que chaque axe du faisceau est perpendiculaire à
l’une des six surfaces du fantôme. Les paramètres à optimiser sont les 96 fluences indépendantes qui
sont initialisées à la même valeur. La dose prescrite pour la cible est toujours fixé à 50Gy.

Les deux plans de coupe définis par les axes et pour le premier et les axes et du repère
fantôme pour le second sont présentés ci-dessous. Sur la Figure 87(a) et (a bis), nous constatons que la
dose est homogène dans le volume cible et l’OAR ce qui est dû à l’initialisation. Ceci est confirmé par
l’histogramme dose-volume initiale, illustré en Figure 89 (a). Au bout de la 5ème itération, illustré par
Figure 80(b) et (b bis), le volume cible est mieux irradié tandis que l’OAR est mieux épargné. Au bout
de la 30ème itération (Figure 80(c) et (c bis)), la dose déposée dans le volume cible se rapproche de la
dose prescrite et la dose continue à diminuer dans l’OAR. Au bout de la 200ème itération, la dose
déposée au bord du volume cible atteint 50Gy et la dose déposée dans l’OAR est quasi nulle. Ceci est
réalisé grâce à la modulation des faisceaux RCMI.

Plan de coup Y - Z à l'initialisation Plan de coup X - Y à l'initialisation

0.06

-4 -4 0.05

-2 -2 0.04
X (cm)
Z (cm)

0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)
(a) (a bis)

Plan de coup Y - Z à 5ème itération Plan de coup X - Y à 5ème itération

0.06

-4 -4 0.05

-2 -2 0.04
X (cm)
Z (cm)

0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)
(b) (b bis)

114
 Plan de coup Y - Z à 30ème itération Plan de coup X - Y à 30ème itération
0.06

-4 -4 0.05

-2 -2 0.04

X (cm)
Z (cm)

0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)
(c) (c bis)

Plan de coup Y - Z à 200ème itération Plan de coup X - Y à 200ème itération

0.06

-4 -4 0.05

-2 -2 0.04
X (cm)
Z (cm)

0 0 0.03

2 2 0.02

4 4 0.01

6 6 0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y (cm) Y (cm)
(d) (d bis)

Figure 87 Les plans de coup: fantôme homogène d’eau avec la cible de forme ‘C’, faisceau type RCMI sans variation
d’incidence

Pour un faisceau de chaque plan de coupe, nous pouvons comparer sa modulation initiale et finale. Sur
la Figure 88(a), la carte de fluence est homogène au départ, elle est composée de 16 fluences
uniformes. Sur la Figure 88(b), la carte s’évolue de manière intelligente, elle a évité d’irradier l’OAR
qui est entouré par le pointillé rouge. Sur la Figure 88(c), le faisceau est modulé de telle sorte qu’il
n’éclaire quasiment que le volume cible. Sur le Figure 88(d), le faisceau s’est conformé à la forme en
‘C’ de la cible en enlevant les faisceaux élémentaires proches du OAR.

Modulation faisceau initial Modulation faisceau N1

0.1

0.09
-4 -4

0.08

-2 0.07 -2

0.06
0 0
X(cm)

0.05

0.04
2 2

0.03

4 0.02 4

0.01
6 6
0
-4 -2 0 2 4 6 -4 -2 0 2 4 6
Z(cm)

(a) (b)

115
 Modulation faisceau N4 Modulation faisceau N5

-4 -4

-2 -2

0 0
X(cm)

Y(cm)
2 2

4 4

6 6

-4 -2 0 2 4 6 -4 -2 0 2 4 6
Y(cm) Z(cm)

(c) (d)

Figure 88 L’évolution de la modulation du faisceau selon différents plans de coup

Histogramme dose-volume Histogramme dose-volume

1 1
Cible Cible
0.9 OAR 0.9 OAR

0.8 0.8

0.7 0.7
Volume [%]

Volume [%]

0.6 0.6

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.01 0.02 0.03 0.04 0.05 0.06
Dose Dose
(a) (b)

Figure 89 Les histogrammes dose-volume initiale (a) et finale (b) : fantôme homogène d’eau avec la cible de forme ‘C’,
faisceau type RCMI sans variation d’incidence

Sur la Figure 89 (b), nous voyons que la totalité du volume cible a reçu une dose supérieure à 33Gy.
En outre, 52% du volume cible a reçu la dose prescrite (50Gy) sachant que cette dose n’est présente
que dans 10% du volume cible à l’initialisation.

Pour l’OAR, la dose déposée est beaucoup diminuée sachant que l’OAR reçoit une forte dose à
l’initialisation.

Si on compare ce plan avec celui de type Cyberknife (résumé par la Figure 85), on peut constater que
le plan de type RCMI atteint mieux l’objectif : sa conformité et son homogénéité dans le volume cible
sont légèrement améliorées sans avoir surexposé l’OAR. Ceci est confirmé aussi par la fonction
objectif, (Figure 90). La fonction converge plus vite (convergence après 40 itérations contre 200
itérations) vers une valeur plus petite que celle du Cyberknife (0.1·10-3 contre 0.25·10-3).

116
 x 10
-3 Fonction objectif
2.5

1.5

0.5
X: 30
X: 200
Y: 0.0001446
Y: 0.000101

0
0 20 40 60 80 100 120 140 160 180 200
Nombre d'iteration

Figure 90 La fonction objectif : fantôme homogène d’eau avec la cible de forme ‘C’, faisceau type RCMI sans
variation d’incidence

3.3.2 Plan isocentre (23 mins)

Dans cette section, nous supposons que le plan de traitement possède un isocentre sur lequel les
axes de chaque faisceau RCMI se croisent. Dans un premier temps, nous considérons que le point
isocentre coïncide avec l’origine du repère fantôme hétérogène. Ceci introduit une contrainte
supplémentaire pour le déplacement du point d’origine du faisceau RCMI puisque alors
, ce qui signifie que dépend totalement de . Les calculs détaillés du gradient est en
annexe E.

La Figure 91 nous montre la comparaison de l’histogramme dose-volume au bout de 30 itérations avec

celui de départ. Nous pouvons dire que l’optimisation a modifié le plan dans le bon sens en tenant
compte de cet isocentre. Ce plan a été légèrement amélioré par rapport au plan précédent qui
n’autorise que la variation des modulations des faisceaux. Ceci est chiffré par la valeur de la fonction
objectif à la 30ème itération, illustré par les Figure 90 et Figure 92 (0.1446·10-3 contre 0.1402·10-3).

Histogramme dose-volume Histogramme dose-volume

1 1
Cible Cible
0.9 OAR 0.9 OAR

0.8 0.8

0.7 0.7 X: 0.048

Y: 0.6158
Volume [%]

Volume [%]

0.6 0.6

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.01 0.02 0.03 0.04 0.05 0.06
Dose Dose
(a) (b)

Figure 91 Les histogrammes dose-volume initial (a) et final (b) : fantôme homogène d’eau avec la cible de forme ‘C’,
faisceau type RCMI Isocentrique

117
 x 10
-3 Fonction objectif
2.5

1.5

0.5
X: 30
Y: 0.0001402

0
0 5 10 15 20 25 30
Nombre d'iteration

Figure 92 La fonction objectif : fantôme homogène d’eau avec la cible de forme ‘C’, faisceau type RCMI Isocentrique

3.3.3 Plan 3D libre (76 mins)

Dans cette section, nous libérons totalement le choix des incidences des faisceaux RCMI. Pour chaque
faisceau RCMI, les 5 paramètres géométriques et 16 fluences sont à optimiser. La dose prescrite pour
la cible est toujours fixée à 50Gy.

L’HDV obtenu au bout de 30 itérations est illustré en Figure 93(b). Nous pouvons constater que ce
plan est moins bon par rapport au plan obtenu dans la section précédente. À titre d’exemple, ce plan a
inclus 60,86% du volume cible dans l’isodose 96% de la dose prescrite au lieu de 61,58% dans le plan
précédent. Ceci est aussi confirmé par une augmentation de la fonction objectif à la 30ème itération
(0.1431·10-3), illustré en Figure 94(a). En revanche, si on pousse le calcul jusqu’à 100 itérations, nous
pouvons constater une légère amélioration tel que 61,92% du volume cible est dans l’isodose 96%.
Ceci est aussi confirmé par la fonction objectif (0.0926·10-3), illustrée en Figure 94. Ainsi, cette valeur
que l’on obtient au bout de 100 itérations est meilleure que celle obtenue avec le plan angle fixe au
bout de 200 itérations (0.101·10-3). Par conséquent, nous pouvons conclure que ce plan ayant le plus
de degré de liberté converge moins vite que le plan isocentre, ceci est dû au nombre important des
paramètres indépendants. Mais le fait qu’il y ait plus de degrés de liberté permet une solution finale
meilleure (à 100 itérations). Donc, en partant sur une balistique définie par un médecin, nous pouvons
faire bénéficier au plan de traitement de ce petit ajustement des incidences des faisceaux pour trouver
une solution optimisée.

Histogramme dose-volume Histogramme dose-volume

1 1
Cible Cible
0.9 OAR 0.9 OAR

0.8 0.8

0.7 0.7 X: 0.048

Y: 0.6086
Volume [%]

Volume [%]

0.6 0.6

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.01 0.02 0.03 0.04 0.05 0.06
Dose Dose
(a) (b)

Figure 93 Les histogrammes dose-volume initiale (a) et 30ème itération (b) : fantôme homogène d’eau avec la cible de
forme ‘C’, faisceau type RCMI avec variation d’incidence

118
 x 10
-3 Fonction objectif Histogramme dose-volume
2.5 1
Cible
0.9 OAR

2 0.8

0.7 X: 0.048
Y: 0.6192

Volume [%]
1.5 0.6

0.5

1 0.4

0.3

0.5 0.2
X: 30
X: 100
Y: 0.0001431 0.1
Y: 9.268e-005

0 0
0 10 20 30 40 50 60 70 80 90 100 0 0.01 0.02 0.03 0.04 0.05 0.06
Nombre d'iteration Dose
(a) (b)

Figure 94 La fonction objectif (a) et histogramme dose-volume finale (b) : fantôme homogène d’eau avec la cible de
forme ‘C’, faisceau type RCMI avec variation d’incidence

4 Fantôme hétérogène Penelope Poumon-Os

4.1 Conditions de simulation
Nous avons validé notre méthode d’optimisation sur un fantôme homogène dans la section précédente,
maintenant, nous allons appliquer cette méthode sur un fantôme hétérogène.

Le fantôme poumons-os est de taille 27×24×64 cm3, il est constitué de 3 volumes simples représentant
deux parallélépipèdes équivalent poumons et une colonne équivalent os représentant la colonne
vertébrale. Le reste du fantôme est rempli d’eau pour modéliser le corps humain. Les trois
composantes d’hétérogénéités sont modélisées par des mailles parallélépipédiques rectangles. Le
fantôme est présenté sur la Figure 95. Chaque poumon coloré en jaune est de taille 6×9×16 cm3,
l’ensemble se situe au centre du fantôme. La colonne vertébrale colorée en gris a une taille de 3×3×16
cm3 La tumeur colorée en rouge de taille 3×2×3 cm3 se situe en bas à droite du poumon gauche.

(a) (b)

Figure 95 Le fantôme hétérogène poumons-os : (a)Vue du plan de coup Y-Z (b) Vue du plan de coup X-Y

Trois faisceaux RCMI de section 4×4 cm2 sont utilisé (16 faisceaux élémentaires par faisceau RCMI).
On autorise les incidences des faisceaux à évoluer indépendamment pendant l’optimisation. En tenant
compte de la modulation des faisceaux, nous avons 63 (3×(5+16)) paramètres à optimiser.

119
Pour définir la fonction de coût, nous allons identifier les différents volumes d’intérêt : la cible et les
OAR. Comme nous l’avons précisé ci-dessus, la tumeur est imbriquée dans le poumon gauche. Nous
allons minimiser l’écart entre la dose prescrite et la dose déposée. Le premier OAR contient les
poumons sains et le second la colonne vertébrale contenant la moelle épinière.

Selon la littérature, la dose de tolérance du poumon dépend fortement, mais pas uniquement, du
volume irradié. En clinique, si le poumon dans son ensemble ne peut recevoir qu’une dose faible, des
petites proportions de l’organe tolèrent des doses élevées. Dans un premier temps, nous avons choisi la
dose pulmonaire moyenne en tant que critère dosimétrique pour constituer la fonction de coût. En
outre, on considère les deux poumons comme un seul organe pair et l’estimation de la dose reçue est
réalisée sur le volume pulmonaire total auquel est soustrait le volume cible.

La dose maximale habituellement admise pour la moelle épinière est de l’ordre de 45–50 Gy
(Habrand, et al., 2010). Dans le but d’établir une fonction de coût facilement dérivable, nous avons
choisi de minimiser la dose déposée dans la colonne vertébrale sans faire intervenir une dose de
tolérance précise.

La fonction de coût est quadratique sous la forme :

Où :
est la dose prescrite à la cible,
est la dose moyenne autorisée par le poumon,
représente la dose déposée dans le ième voxel du volume cible,
représente la dose déposée dans le jème voxel du poumon sain,
représente la dose déposée dans le kème voxel de la colonne vertébrale,
, et représentent respectivement les nombre de voxels du volume cible, du poumon sain et de
la colonne vertébrale.

4.2 Échantillonnage des points de contrôles

Avec une résolution de 2 mm, le nombre de voxels dans le fantôme entier s’élève à environ 5 millions.
Même restreint aux volumes d’intérêt (Cible et OAR), ce nombre est toujours significatif : 216 000
voxels. Dans le cadre de notre étude, nous avons réduit ce nombre en définissant des points de
contrôle pour accélérer le calcul. Cela revient à du sous-échantillonnage des voxels.

Les techniques d’échantillonnage pour les critères dosimétriques ont déjà été étudiées par (C Martin, et
al., 2007). L’idée de base est de sélectionner une partie des voxels de chaque volume d’intérêt, ces
voxels sélectionnés – et uniquement eux – vont contribuer au calcul de la fonction objectif. La
sélection peut se faire avec un taux fixe pour chaque volume ou un taux variable selon l’importance de
ce dernier. La loi de probabilité utilisée pour la sélection est la loi uniforme.

En pratique, nous avons choisi une technique moins sophistiquée. Pour les deux poumons, on
sélectionne un voxel sur deux selon les 3 dimensions, ce qui en nombre de voxels considérés est
équivalent à un changement de résolution de 2 mm à 4 mm. Pour la colonne vertébrale, nous ne nous
intéressons qu’à la partie proche de la cible (la moitié de la colonne vertébrale) avec une résolution de
4 mm également. Par conséquent, le taux d’échantillonnage pour les poumons est 12.5% et 6.25%
pour la colonne vertébrale.

4.3 Étude de cas (82 mins)

120
Dans cette section, nous allons analyser le résultat d’optimisation dans le fantôme hétérogène. La dose
prescrite à la cible est 80 Gy (correspond à 0.08 sur l’échelle de couleur). Un pas fixe est utilisé.

Sur la Figure 96, nous affichons les plans de coupe (définis par les axes et ) à différentes itérations
et leur HDV correspondant. Le volume cible est entouré par un rectangle en pointillé rouge, l’os est
entouré par un carré en pointillé blanc à gauche, les deux poumons sont entourés par deux rectangles
en pointillé blanc. La dose n’est calculée que dans ces volumes d’intérêt.

Nous pouvons constater que, à l’initialisation, l’intersection des trois faisceaux se situe dans le
poumon du haut mais en dehors de la cible : la cible ne reçoit pas une dose suffisante et le poumon en
haut est trop irradié (au moins pour ce plan de coupe). Au bout de 30 itérations, le volume de la cible
recevant une dose suffisante augmente et l’os est mieux épargné, confirmé par l’HDV de la Figure 96
(b bis). Au bout de 60 itérations, l’intersection des trois faisceaux est bien décalée vers la gauche pour
cibler la tumeur, et le volume d’os recevant une faible dose (ex. 10 Gy) a diminué. Au bout de la 90ème
itération, la dose minimale déposée dans la cible est 40 Gy et la moitié du volume reçoit une dose
supérieure à 70 Gy. Le poumon est mieux épargné par rapport à la 60ème itération. Par exemple,
12.02% de volume (Figure 96 (d bis)) reçoit une dose supérieur à 10 Gy contre 12.74% (Figure 96 (c
bis)) à la 60ème itération. Sachant que les parties saines des poumons représentent 30 400 voxels, cette
différence en pourcentage représente 218 voxels pour notre fantôme. La dose est passée sous le seuil
pour ces 218 voxels.
Plan de coup X - Y à l'initialisation Histogramme dose-volume
0.08 1
-12 Cible
0.9 Poumon
0.07
-9 Os
0.8
0.06
-6
0.7
0.05
Volume [%]

-3
0.6
X (cm)

0 0.04 0.5

3 0.4
0.03

6 0.3
0.02
0.2
9
0.01
0.1
12

0 0
-9 -6 -3 0 3 6 9 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
Y (cm)
Dose
(a) (a bis)

Plan de coup X -Y à 30ème itération Histogramme dose-volume

0.08 1
-12 Cible
0.9 Poumon
0.07
-9 Os
0.8
0.06
-6
0.7
0.05
Volume [%]

-3
0.6
X (cm)

0 0.04 0.5

3 0.4
0.03

6 0.3
0.02
0.2
9
0.01
0.1
12

0 0
-9 -6 -3 0 3 6 9 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
Y (cm)
Dose
(b) (b bis)

121
 Plan de coup X - Y à 60ème itération Histogramme dose-volume
0.08 1
-12 Cible
0.9 Poumon
0.07
-9 Os
0.8
0.06
-6
0.7

Volume [%]
-3 0.05
0.6
X (cm)

0 0.04 0.5

3 0.4
0.03

6 0.3
0.02
X: 0.01
0.2
9 Y: 0.1274
0.01
0.1
12

0 0
-9 -6 -3 0 3 6 9 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
Y (cm)
Dose
(c) (c bis)

Plan de coup X - Y à 90ème itération Histogramme dose-volume

0.08 1
-12 Cible
0.9 Poumon
0.07
-9 Os
0.8
0.06
-6
0.7

-3 0.05 Volume [%] 0.6

X (cm)

0 0.04 0.5

3 0.4
0.03

6 0.3
0.02
0.2 X: 0.01
9 Y: 0.1202
0.01
0.1
12

0 0
-9 -6 -3 0 3 6 9 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
Y (cm)
Dose
(c) (d bis)

Figure 96 Les plans de coup et les HDVs : fantôme hétérogène, faisceau type RCMI

Figure 97, nous pouvons voir la décroissance de la fonction objectif, la pente de la décroissance
diminue au fur à mesure des itérations. Néanmoins, la convergence n’est pas encore atteinte au bout de
90 itérations. Dans le cadre de cette étude, un pas fixe unique est défini pour tous les paramètres
angulaires. Nous avons constaté que, pendant l’optimisation, la dynamique du gradient pour chaque
paramètre n’est pas prise en compte. Ceci ralenti la convergence de l’algorithme. Dans le but
d’améliorer la performance de l’algorithme, nous pouvons envisager de choisir un pas distinct et
variable pour chaque paramètre en fonction de leur dynamique au court de chaque itération.

x 10
-3 Fonction objectif
7

6.5

5.5

4.5

3.5

3
0 10 20 30 40 50 60 70 80 90
Nombre d'iteration

122
 Figure 97 La fonction objectif : fantôme hétérogène, faisceau type RCMI

D’ailleurs, plus généralement, les observations de plusieurs autres simulations faites mais non
présentées ici car n’apportant rien de probant en plus, nous ont permis de constater que la méthode de
descente de gradient tombe facilement sur des minimums locaux. L’état de l’art indique que ce
problème peut être amélioré par :

 Une meilleure initialisation du plan : ceci est cohérent avec le principe initial qui vise une
optimisation d’un plan de traitement proposé par le physicien.
 Un échantillonnage plus sophistiqué (ex. échantillonnage stochastique où les voxels
considérés d’une itération à l’autre ne sont pas les mêmes).
 La méthode d’optimisation du gradient conjugué

Le temps d’optimisation est de 3.7 minutes pour le cas de RCMI avec une balistique fixe. Ceci est
logique car la propagation du faisceau dans le fantôme a fait une fois pour tous, la dose calculé par
chaque faisceau élémentaire d’intensité unitaire est constante. Ceci implique l’augmentation
importante du temps de calcul pour les autres cas de test qui nécessitent un recalcul complet de la dose
(propagation puis évaluation de la dose car nouvelle balistique) à chaque itération. Le temps de
recalcul dépend de la complexité du maillage dans le fantôme (le maillage influence le nombre de
SBIM à propager, mais aussi du coup le nombre de SBIM participant à la somme pondérée pour le
calcul total de dose). Ceci explique le fait que le temps de calcul pour le fantôme homogène d’eau est
plus faible que celui du fantôme hétérogène. Plus précisément, le cas du fantôme hétérogène nécessite
82 minutes de calcul pour 90 itérations, chaque itération faisant le calcul de dose de 48 faisceaux fins
(3 faisceaux RCMI de 16 faisceaux élémentaires) et autant de calcul de gradients, dans un fantôme
composé de 288 mailles. Par rapport au 20 minutes ciblées, les conditions et sources d’évolutions des
temps de calcul par rapport à ce code nous laissent penser que Doséclair peut effectivement atteindre
cet objectif.. Le code s’exécute en mono-cœur avec la boucle d’optimisation implémentée sous
MATLAB et appelant les fonctions de calcul Doséclair comme mexfiles. Ces fonctions étaient
implémentées sur une structure de classes initiée en 2006 pour être particulièrement modulaire afin de
tester toutes sortes de modèles (via principalement des mécanismes d’héritage), au détriment
naturellement des performances qu’un modèle unique pourrait atteindre en étant codé de façon
optimale.

Le temps de calcul de 82 minutes pour une telle configuration est donc finalement raisonnable. Un
simple passage en hexa ou octo-cœurs ferait déjà passer la simulation en dessous des 20 minutes. Si on
fait des projections sur l’augmentation des calculs sur des fantômes plus réalistes (donc avec plus de
mailles), on peut vite imaginer avoir un rapport 100 à 1000 mais il s’avère qu’empiriquement ce sont
aussi les ordres de grandeurs que l’on peut gagner en passant de ces codes génériques en C++ et
MATLAB à une implémentation optimisée de l’algorithmes faites par un programmeur expérimenté.
Le bilan serait alors neutre et l’objectif de temps atteint. Dans le cadre d’une preuve de concept avant
justement d’approfondir les développements, les résultats de cette thèse nous semblent plutôt positifs.

À ce jour, nous avons vu que la méthode d’optimisation fondée sur la technique de calcul de dose
rapide Doséclair nous permet, dans un temps de calcul raisonnable, d’améliorer le plan de traitement
de façon continue en les paramètres, et en 3 dimensions. Ces résultats obtenus ouvrent de nombreuses
perspectives dans le domaine de l’optimisation des plans de traitement en radiothérapie.

123
Conclusion et perspectives
L’objectif de cette thèse était de montrer, sur des cas simples de validation, la faisabilité d’application
du nouvel algorithme de calcul de dose et son utilisation pour l’optimisation du plan de traitement.

La première phase du travail a consisté la prise en main de la technique de calcul de dose et du code de
calcul Doséclair. Une reformulation tensorielle de la propagation du faisceau dans le fantôme a été
effectuée. Les paramètres du faisceau sont redéfinis par les angles d’Euler. Les codes sources de
Doséclair ont été modifiés pour prendre en compte ces nouveaux paramètres et cette formulation
tensorielle pendant la propagation du faisceau. La modification la plus importante est le module
‘simulation’ du code de calcul.

Dans le but de valider ce calcul tensoriel dans les configurations qui sont présentés lors les travaux de
Baptiste Blanpain, quatre tests ont été réalisés. Les distributions de dose obtenues avec la formulation
tensorielle sont identiques aux anciennes.

Cette reformulation ouvre la possibilité de définir le gradient de la dose en fonction des nouveaux
paramètres de la balistique de traitement. Une formulation tensorielle du gradient de dose et sa
propagation dans le fantôme a été effectuée. Cette formulation a été ensuite implémentée dans
Doséclair.

Trois types de test ont été présentés pour la validation des calculs de gradients : test individuel, test du
fantôme homogène, test du fantôme hétérogène. La méthode de descente de gradient avec une fonction
objectif quadratique classique est utilisée pour cette validation. Les HDVs et les plans de coupe ont été
utilisés pour évaluer le plan final. La convergence des fonctions objectifs a été observée dans la
majorité des simulations.

L’objectif initial de la thèse a donc été atteint au sens où la preuve de concept est faite et qu’aucun
résultat ne montre de caractéristiques rédhibitoires (temps de calculs, précision). Toutefois la méthode
d’optimisation développée n’a pour l’instant pu être testée que sur des cas simples. Il reste encore
beaucoup de travaux avant de rendre cette méthode opérationnelle sur les cas cliniques :

La segmentation du fantôme en mailles homogènes sur des données réelles sera à évaluer pour trouver
le compromis entre la taille (et donc le nombre) des mailles de forme parallélépipède rectangle et la
précision de la quantification que cela implique sur les organes réels. Une étude sur la segmentation
avec des surfaces quadratiques sera intéressante. Dans un premier temps, les géométries quadratiques
utilisées dans les codes de Monte Carlo peuvent servir comme composants de références. Le
découpage du faisceau en SBIM et la propagation du faisceau pourront y être testés par une approche
similaire à celle suivie dans la thèse.

Le réglage des fonctions de pondération est basé sur des données de simulation Monte Carlo. La forme
de la fonction est choisie en fonction de sa capacité de modéliser le dépôt de dose dans les zones de
déséquilibre électronique. La dose diffusée est géré par le scaling latéral. Ce coefficient de scaling est
également obtenu expérimentalement avec la dose de référence Monte Carlo. Très certainement, ces
modèles devront être remis en cause et revalidés voire modifiés lors du passage aux données réelles.
Toutefois le principe ne devrait pas changer sur le fond, puisque l’on peut en donner des explications
physiques et une validation numérique par comparaison avec les codes de Monte-Carlo qui sont les
codes de calculs les plus précis. La même question se pose pour le modèle de spline utilisé pour les
modèles homogènes, la précision devant être validée en cas clinique. Toutefois là encore, la bonne
précision avec les simulations de Monte-Carlo est plutôt prometteuse.

La méthode de descente de gradient à pas fixe est utilisée pour prouver que le plan de traitement
pourra être optimisé à l’aide des gradients de dose calculé au cours de la propagation. Les résultats

124
nous montrent que les plans sont améliorés au fur à mesure des itérations. Sur des cas réels, cette
méthode qui est la plus simple de toutes les descentes de gradient montrera surement ses limitations
classiques : convergence lente, minimums locaux nombreux. D’autres algorithmes de descentes de
gradient plus sophistiqué seront à évaluer, comme le gradient conjugué. Parmi les pistes on peut même
imaginer, tout comme on a pu dériver une fois le calcul de dose pour calculer le gradient, dériver
partiellement une seconde fois l’ensemble et avoir une méthode de calcul du hessien par propagation.
S’ouvre alors même des méthodes comme Newton, le pas optimal ou Levenberg-Marquardt.

Pour accélérer le temps de calcul, un échantillonnage est réalisé dans le test fantôme hétérogène. Une
étude précise sur l’influence de l’échantillonnage (nombre de point, densité d’échantillon par zone…)
sur la contribution de la fonction de coût sera intéressante. Une politique d’échantillonnage efficace
devra être appliquée dans un traitement réel. Des études théoriques sont abordées par plusieurs
auteurs. Un large d’étude sur ce sujet est disponible.

Le calcul de gradient analytique ouvre la voie à l’optimisation de trajectoire dans le cas du Cyberknife
et plus généralement de l’arcthérapie. Pour cela, on pourrait considérer en première approximation que
l’intégration de la dose sur la trajectoire serait égale à la somme des doses avec des points d’origine
régulièrement répartie le long de la trajectoire. En supposant que cette trajectoire soit définie par un
petit nombre de valeurs références (d’une équation) ou de points références (les brisures sur une
trajectoire en ligne brisée), on peut aisément établir les relations entre ces références et tous les
paramètres balistiques desdits faisceaux intermédiaires. Un gradient reliant la dose déposée sur la
trajectoire à ces références est donc calculable ou pourrait être exploité pour optimiser tous les
paramètres de la trajectoire, avec un coût en puissance de calcul équivalente à celle pour le RCMI.

La prise en compte d’un faisceau poly-énergétique utilisé en clinique autre que le faisceau mono-
énergétique devrait se faire naturellement, en faisant apprendre les modèles homogènes par un
faisceau poly-énergétique. Le durcissement du faisceau serait pris en compte naturellement puisqu’il a
lieu aussi en milieu homogène et que le modèle de propagation reposent sur une mise à l’échelle entre
matériaux n’est pas affecté par l’énergie du faisceau. On évalue même que les incertitudes,
correspondant à l’intégration des incertitudes sur le durcissement du spectre fois les légers écarts de
dépôt de dose selon l’énergie, seront largement inférieures aux erreurs dues aux incertitudes initiales
sur le spectre lui-même. Ceci restera à être vérifié, et la plateforme Doséo sera un parfait outil pour
cela.

Au-delà du domaine d’application de la radiothérapie, l’esprit d’un modèle analytique à l’échelle

d’une maille et sa propagation automatique pourraient aussi être intéressant dans d’autre domaine
d’application, par exemple le calcul de transfert de l’énergie dans le bâtiment, le calcul de transfert de
chaleur dans les réacteurs nucléaires…

125
Annexe
Tableau de notation :

Scalaire …
Constante …
Fonctions
Point de l’espace …
Vecteur …
Matrice …
Tenseur d’ordre 2

Coefficients d’un tenseur d’ordre 2

Tenseur d’ordre 4
Coefficients d’un tenseur d’ordre 4
Numéro du SBIM en exposant dans la
propagation axiale
Matériaux du milieu en exposant dans la
propagation latérale
Distance entre deux points géométriques
Norme d’un vecteur
Produit scalaire entre deux vecteurs par exemple
Indice faisceau élémentaire d’un faisceau
RCMI
Indice faisceau élémentaire d’un ensemble
de faisceau RCMI
Repère
Coordonnées en repère fantôme (repère de
référence ) du point

Et du vecteur

Coordonnées en repère faisceau ( ) du

point

Coordonnées en repère homogène ( ) du

point

Tenseur d’ordre 1 correspondant aux

coordonnées d’un point étendues d’un 1 afin
de pouvoir appliquer les formules ,
tensorielles
Gradient de dans la direction

Dérivé de par rapport à

126
Annexe A Propagation axiale du tenseur de projection dans les mailles support du fantôme

Nous cherchons à caractériser le tenseur de passage de à , noté , et défini par :

Soit un point ayant les coordonnées dans le repère fantôme , ses coordonnées dans le
repère faisceau peuvent être déduites par la relation tensorielle :

Dans la deuxième partie, Chapitre IV, section 1, nous avons déjà montré que:

La composante (resp. ) est le produit scalaire de l’axe (resp. ) et le point d’origine du

faisceau. Au final, seule la composante varie d’un SBIM à l’autre. Nous allons voir comment
calculer à partir de .

Figure 98 Propagation axial.

Nous avons vu dans la deuxième partie, Chapitre IV section 2.1 que la profondeur radiologique du
SBIM est calculée de manière récursive :

(A.1)

Pour un point dans le SBIM , la profondeur projetée dans le repère faisceau, , peut
s’écrire :
(A.2)

127
Pour calculer la profondeur projetée dans le repère homogène, , on doit tenir compte de l’offset
, la distance parcourue depuis le point d’origine du faisceau selon la direction pour rentrer
dans le SBIM , et de la profondeur radiologique du SBIM , notée .

D’abord, nous calculons la profondeur du point dans son SBIM , qui vaut .
Ensuite, nous ajoutons, à cette profondeur, la profondeur radiologique du SBIM pour avoir la
profondeur projetée du point dans le repère homogène, d’où l’équation suivante :

(A.3)

En remplaçant de l’équation (A.3) par l’équation (A.2), on obtient :

(A.4)

Or

Nous pouvons en déduire que :

(A.5)

En combinant (A.1) et (A.5), on obtient :

(A.6)

Or ,

L’équation (A.6) devient :

(A.7)

Après une factorisation de l’équation (A.7), cette dernière nous donne la formule de propagation de la
composante du SBIM au SBIM :

+ + (A.8)

Notons que pour le SBIM , , qui peut être également déduit de (A.8)
en posant .

Une fois que nous savons comment calculer à partir de , nous pouvons déduire de
l’ensemble des formules le tenseur de propagation de projection d’ordre 4, ( en indice
désigne projection), pour passer de à . Ce passage revient à un produit tensoriel doublement
contracté :

Comme la grande majorité des composantes sont constantes, le tenseur de propagation est basé sur le
tenseur identité. Seul le nouveau coefficient requière plus de calculs. Pour cela il faut voir
comment calculer l’offset . Rappelons-nous que est la distance parcourue depuis le point

128
d’origine du faisceau selon la direction pour rentrer dans le SBIM par n’importe quelle
surface. Etant donnée le point d’origine du faisceau et la direction du faisceau incident , pour un
SBIM , on cherche à connaitre par quelle surface rentre le faisceau et la distance parcourue par ce
dernier.

D’abord, quelques notations pour la suite des calculs :

(resp. , ) : Points d’intersection du faisceau et l’une des deux surfaces perpendiculaires à

l’axe (resp. , ) de la maille support du SBIM . Le point d’entrée du faisceau dans le SBIM ,
, appartient à l’ensemble de , et .

resp. , ): Distance parcourue depuis le point d’origine du faisceau selon la direction

pour traverser la plus proche des deux surfaces perpendiculaires à l’axe (resp. , ) de la maille
support du SBIM , appartient à l’ensemble de , et .

, (resp. , , , ): L’abscisses (resp. l’ordonnées, l’altitudes) des deux

surfaces perpendiculaires à l’axe (resp. , ) de la maille support du SBIM .

Figure 99 Calcul de λ<s>

Sur l’exemple de la Figure 99, le point est l’intersection de l’axe du faisceau et la surface en haut
de la maille support du SBIM .

Par définition, le point est sur l’une des deux surfaces perpendiculaires à l’axe, son altitude dans
le repère fantôme, est l’une des deux valeurs et , ce qui implique qu’elle est bornée par
ces deux valeurs :

≤ ≤ (A.9)

Le vecteur peut s’exprimer en fonction du :

= (A.10)

Selon l’axe du fantôme, l’égalité (A.10) peut encore s’écrire comme:

(A.11)

En remplaçant de l’inégalité (A.9) par l’équation (A.11), on obtient :

≤ ≤ (A.12)

129
En fonction du signe de , l’inégalité (A.12) peut être transformée :

En suivant le même raisonnement, on peut en déduire les deux relations suivantes pour les autres
axes :

Selon les diverses configurations, la distance parcourue pour entrer dans la maille support du SBIM
est représentée par les bornes inférieures . La distance parcourue pour sortir de la
maille support du SBIM est représentée par les bornes supérieures .

Géométriquement, la plus grande valeur positive ou nulle, parmi les trois bornes inférieures,
représente la vraie distance d’entrée dans la maille. La plus petite valeur, parmi les trois bornes
supérieures, représente la vraie distance de sortie. D’où

= (A.13)

L’équation (A.13) nous permet de calculer la distance d’entrée dans la maille support du SBIM en
combinant les fonctions min et max.

Figure 100 Découpage d’un faisceau en deux SBIMs

Dans la réalité, il arrive que le faisceau rencontre des interfaces obliques, que ce soit à son entrée dans
le fantôme, ou en profondeur. Quand l’angle du faisceau par rapport à ces interfaces est grand, nous
avons proposé de séparer le faisceau incident latéralement, en plusieurs SBIMs, illustré par la Figure
100.

Dans ce cas, nous avons besoin de et pour la création des deux tenseurs de projection. En
pratique, le barycentre de la section d’entrée du SBIM , notée , est propagé pendant la
propagation du faisceau dans le fantôme. Pour chaque SBIM , ce centre est projeté sur le plan de la
section du faisceau incident selon la direction pour calculer le point d’origine qui est
ensuite utilisé pour calculer le associé.

130
Annexe B Propagation latérale du tenseur de projection dans les mailles voisines du
fantôme

Figure 101 Propagation latérale du tenseur de projection dans un fantôme hétérogène

Sur la Figure 101, le tenseur de projection, , est calculé pour le SBIM de la maille d’eau. Ce
tenseur est propagé de la maille d’eau à sa maille voisine d’os en tenant compte du changement de
matériau.

Prenons l’exemple du point , nous cherchons à connaitre sa projection dans le repère

homogène d’eau grâce au nouveau tenseur de la maille d’os, selon = .

Dans la propagation latérale, on gère seulement la diffusion latérale, ce qui implique que la profondeur
projetée du point par le tenseur et par le tenseur est identique dans le repère
homogène. Cette profondeur projetée dans le repère d’homogène d’eau, est calculée grâce
à la troisième ligne, indexé par , du tenseur de projection , expliqué en Annexe A. On note :

= =

Nous nous focalisons donc sur le calcul des deux composantes latérales projetées, et
du point dans le modèle homogène d’eau.

On note, , un vecteur géométrique dans le fantôme hétérogène, et , ce vecteur

projeté dans le modèle homogène du milieu en tenant compte de l’hétérogénéité.

On projette le point orthogonalement sur l’axe du faisceau, est le point projeté, est
l’intersection de cette projection et l’interface des deux mailles. Donc, est perpendiculaire à
.

Le vecteur projeté dans le modèle homogène l’eau en tenant compte de l’hétérogénéité,

, peut s’écrire :

= + (B.1)

Or le vecteur géométrique est intégralement dans la maille d’eau, donc sa projection,

, est identique. Le vecteur géométrique est intégralement dans la maille d’os, on
suppose que son projeté lui est proportionnel selon un coefficient :

= (B.2)

est le coefficient de scaling qui dépend des deux matériaux, = / , pour prendre en
compte ce changement de milieu. L’équation (B.1) devient :

131
 = + (B.3)

Or = - , (B.3) peut aussi être transformée :

= - + = +

Finalement, nous obtenons :

= +

On projette sur l’axe du repère d’homogène d’eau :

(B.4)

est la composante latérale du vecteur projeté dans l’eau selon l’axe , noté
. Elle est obtenue grâce à la première ligne, indexée par x, du tenseur de projection de la
maille voisine d’os . Donc, = .

est la composante latérale géométrique du vecteur selon l’axe sans tenir

compte de l’hétérogénéité. Elle est obtenue grâce à la première ligne du tenseur de projection de la
maille support d’eau . Cette composante est noté .

L’équation (B.4) se réécrit :

= (B.5)

Avec le même raisonnement, si on projette sur l’axe du repère d’homogène d’eau, ce

qui correspond aux calculs avec la deuxième ligne, indexé par y, des deux tenseurs de projections ,
l’équation (B.5) devient :

= (B.5 bis)

Nous retrouvons des formules linéaires simples qui permettront de définir facilement un tenseur
additif de propagation latéral, moyennant d’exprimer en fonction des paramètres connus
du fantôme hétérogène et du faisceau incident.

Dans un premier temps, supposons que est dans le plan défini par l’origine du faisceau, l’axe
du faisceau et l’axe du fantôme hétérogène, illustré par la Figure 101. Dans ce plan, on définit un
vecteur unitaire perpendiculaire à . On peut déduire les relations suivantes :

, et sont deux scalaires.

Après simplification,

+ pour (B.6)

132
Or est perpendiculaire à , il est donc porté par ce vecteur unitaire :

= , c est un scalaire (B.7)

On projette sur l’axe du fantôme hétérogène :

, est l’abscisse de l’interface (B.8)

En tenant compte (B.7), l’équation (B.8) peut s’écrire :

On peut déduire :

Donc peut s’écrire :

Par conséquence, peut s’exprimé :

(B.9)

Dans l’équation (B.9), on remplace le vecteur par son expression explicitée en (B.6). Après
simplification, on obtient :

(B.10)

L’expression (B.10) est un produit scalaire d’un vecteur de 4 éléments avec le tenseur d’ordre 1
du point dans le repère du fantôme:

(B.11)

Avec et

On reprend l’équation (B.5) :

On suppose que , est le tenseur de propagation latérale, il est

d’ordre 2. L’équation (B.11) nous permet d’identifier la première ligne du tenseur de propagation
latérale :

133
 ème ème
représente la composante de la ligne et la colonne du tenseur .

En suivant le même raisonnement, la projection de sur l’axe du repère d’homogène

d’eau, , prend la même forme que (B.11), il faut juste remplacer de l’équation
(B.11) par . On obtient le produit scalaire suivant :

Avec et

Il nous permet, avec l’équation (B.5 bis) d’identifier la deuxième ligne du tenseur de propagation
latérale :

Maintenant, nous allons étudier le cas d’une propagation latérale multiple, illustrée par la Figure 102.

Figure 102 Propagation latérale multiples

Supposons que le tenseur de projection de la maille d’eau, , est connu. Ce tenseur est propagé
de la maille d’eau à sa maille voisine d’os par le tenseur de propagation latérale expliqué
précédemment. Nous allons voir comment propager ce tenseur, , à la maille de poumon.

Prenons exemple du point M4, le vecteur projeté dans l’eau, , grâce au tenseur de
la maille de poumon peut s’écrire :

= +

On introduit le tenseur de la maille de poumon :

Et
(B.12)

Si cette projection est réalisée grâce au tenseur de projection de la maille d’os :

= +

134
On introduit le tenseur de la maille d’os :

Et
(B.13)

Or M3 est sur l’interface des deux mailles d’os et poumon, donc (resp.
) et (resp. ) sont identiques. En tenant compte de
(B.13), l’équation (B.12) peut s’écrire :

Et

On peut calculer le tenseur de projection grâce au tenseur de projection , qui est lui-
même calculé à partir du tenseur de projection . La propagation latérale de maille à maille est
donc toujours réalisée par une addition tensorielle, mais le coefficient de scaling dépend des trois
matériaux : celui de la maille support, celui de la maille précédente et celui de la maille courante, selon
.

Dans le cas de la Figure 101, la maille courante est de l’os, la maille précédente et la maille support
sont de l’eau, et on retrouve alors le coefficient .
Cette propagation est exacte pour les points qui sont dans le plan défini par l’axe du fantôme
hétérogène, l’origine et l’axe du faisceau. En pratique, par une approximation, on considère
qu’un point en dehors de ce plan subit le même scaling que sa projection orthogonale dans ce plan.

Nous avons identifié le tenseur de propagation latérale selon la direction du fantôme hétérogène. Par
symétrie, on peut facilement déduire le tenseur de propagation, , selon la
direction et sous la formule générale suivante :

Pour la direction :

, est l’ordonnée de l’interface

135
 Pour la direction :

, est l’altitude de l’interface

En pratique, pour un faisceau incident oblique, différentes zones sont définies dans la maille voisine.
Chaque zone est associée à un tenseur de projection latérale avec une direction de propagation
distincte.

Figure 103 Sélection des surfaces de propagation

Sur l’exemple de la Figure 103, l’axe du faisceau oblique est présenté par le vecteur . Un SBIM
en pointillé rouge est créé à l’entrée du fantôme. Ce SBIM est propagé latéralement dans la maille
de poumon. Pour chaque point de la maille poumon, un scaling sur sa projection latérale est effectué
pour tenir compte de ce changement de milieu.

En pratique, une liste de surfaces candidates est associée au tenseur de propagation,

au moment de la propagation latérale du SBIM. Cette liste est créée par une évaluation des six surfaces
de la maille de poumon. Sur la Figure 103, l’interface commune de la maille d’eau et la maille de
poumon est la surface 1, on note, , le vecteur normal de cette surface. La formule d’évaluation des
surfaces de la maille est :

,est le vecteur normal des 4 surfaces adjacentes de la maille, la sixième

surface orthogonales à étant exclue d’office. Les surfaces vérifiant cette condition sont
sélectionnées.

Dans ce cas, les surfaces 3 ( ) et 4 (>0) sont exclues et la liste des surfaces candidates contient
les surfaces 1 et 2.

136
Ensuite, pour chaque point de la maille, on calcule les normes de avec la
surface 1 et 2. Pour ce point, on compare ces normes et choisi la surface qui donne la norme minimale.
Cette comparaison des normes est équivalente à dire qu’en partant du point et en projetant
orthogonalement ce point sur l’axe , on cherche la surface par laquelle traverse cette projection.
Avec cette surface, nous pouvons calculer définitivement le tenseur de propagation pour le point.

Nous pouvons voir sur la Figure 103, cette sélection de surface résulte d’une partition de la maille de
poumon en deux zones qui sont séparées par la surface encadrée en pointillés noir. Pour les points de
la zone 1 (resp. zone 2) qui présente le volume en haut (resp. en bas), surface numéro 2 (resp. numéro
1) est utilisée pour calculer le tenseur de propagation.

137
Annexe C Gradient du tenseur de projection axial et sa propagation

On cherche à calculer le gradient du tenseur de projection, , par rapport à chaque paramètre

géométrique/balistique . Ce gradient prend lui aussi la forme d’un tenseur de même dimension :

Rappelons-nous que le tenseur de projection axial est sous la forme générale :

= avec ,

La composante (resp. ) est le produit scalaire de l’axe (resp. ) avec le vecteur , où

est le point d’origine du faisceau indexé par . est la composante calculée de proche à
proche pour intégrer la profondeur radiologique. Parmi ces 12 composantes, il n’y a que qui
dépende de la propagation du faisceau dans le fantôme. Les 11 composantes restantes dépendent
seulement des axes (orientation) et du point d’origine du faisceau. Par conséquent, le gradient de ces
11 composantes ne change pas au cours de la propagation.

Pour un faisceau élémentaire indexé par du RCMI indexé par , il est paramétré par =( ;
; ; ; ). Afin de simplifier les écritures, nous ne mettons pas l’indice sous les paramètres
angulaires.

Nous pouvons exprimer le tenseur de projection du faisceau élémentaire :

Or le point d’origine du faisceau indexé par vaut :

Nous pouvons exprimer la composante :

=
=

, origine du faisceau RCMI, est défini par :

138
Et l’axe du faisceau:

Nous pouvons donc exprimer la composante en fonction des cinq paramètres (θ ; ϕ ; α ; β ; γ) :

=

Avec le même déroulement de calcul, nous pouvons exprimer la composante :

=
=
–
( )

Maintenant, nous allons exprimer la composante . Dans l’annexe A, nous avons vu que la
composante peut être calculé de proche à proche par le formule suivant :

+ +

Rappelons-nous que l’exposant désigne le SBIM courant , et désigne le SBIM qui « le suit
géométriquement ». Cette formule nous permet de connaitre le gradient de à partir de celui de
. Plus précisément, nous pouvons écrire :

(C.1)

Où désigne le gradient de ( est un scalaire ou un vecteur) par rapport au paramètre

Nous souhaitons donc pouvoir évaluer le gradient de ce tenseur. Pour presque chaque composante,
nous disposons d’une formule analytique en fonction des paramètres. Nous pouvons alors dériver cette
formule analytique et en déduire également une formule analytique pour le gradient. Cela signifie que
même si le gradient n’est pas le même en chaque point de la maille, on dispose des formules pour
calculer numériquement les composantes de ce gradient en connaissant le point et les paramètres. Ceci
est vrai pour toutes les composantes sauf , pour laquelle nous avions du élaborer notre approche
par propagation. On va donc reproduire ici la même approche.

139
On définit le tenseur de propagation de gradient d’ordre 4, (l’exposant désigne le
paramètre qui peut être n’importe lequel des 5 paramètres ballistique, l’indice désigne gradient),
pour passer de à . Ce passage est réalisé par un produit tensoriel doublement contracté :

Le est composé de 256 (44) composants indexés par deux paires d’indices et ,
, les indices varient de 1 à 4. Ce produit tensoriel doublement contracté est exprimé explicitement:

Nous avons vu que, pour une propagation axiale du tenseur de gradient, est la seul composante qui
change au cours de la propagation, donc, seul les sont non nulles. Or le terme de
l’équation (C.1) peut être identifié comme :

+ +

Nous pouvons déduire, pour un faisceau élémentaire , le tenseur de propagation de gradient

:

L’expression peut s’écrire explicitement :

(C.2)

où désigne

Selon le paramètre , les cinq valeurs numériques calculées par l’équation (C.2) sont constantes au
cours de la propagation. Par conséquent, elles sont précalculées une fois pour toute. Nous allons
expliciter ces cinq valeurs selon :
 , (C.2) devient :

140
  , (C.2) devient :

 , (C.2) devient :

 , (C.2) devient :

 , (C.2) devient :

Dans l’annexe A, nous avons vu que, pour un faisceau élémentaire , peut être calculé par le
formule suivante :

avec

min ,

Nous pouvons déduire que est :

Donc, il faut d’abord évaluer les fonctions min et max pour chaque propagation axiale, ensuite,
calculer le gradient du terme approprié. En pratique, les 6 gradients candidates sont calculés une fois
pour tout. Pendant la propagation, on va simplement choisir, selon le résultat de l’évaluation des
fonctions min et max, quel gradient doit être utilisé pour s’intégrer dans son tenseur de propagation.
Nous allons lister ces six gradients candidats selon :

141
 

On remplace par dans les six expressions précédents

Même principe, on remplace par cette fois ci

Comme les trois vecteurs ne dépendent pas de , les six valeurs sont réduit à trois

Même principe, on remplace par

Nous avons vu jusqu’à présent comment propager de proche à proche le tenseur de gradient grâce au
tenseur de propagation de gradient. Maintenant, nous allons expliciter les cinq tenseurs de gradient de
projection pour les SBIM correspondant à l’entrée du faisceau élémentaire dans le fantôme, noté de
façon générique <1> :

 Tenseur de gradient par rapport à

+
( )

 Tenseur de gradient par rapport à

142
 –

 Tenseur de gradient par rapport à

 Tenseur de gradient par rapport à

 Tenseur de gradient par rapport à

143
Nous avons tous les éléments nécessaires pour initialiser et propager axialement le tenseur de gradient.
Nous allons voir dans l’annexe suivante comment propager latéralement un tenseur de gradient.

144
Annexe D Gradient du tenseur de projection latérale et sa propagation

Le tenseur de projection est propagé latéralement dans la maille voisine par l’addition d’un tenseur de
propagation latérale :

On note le gradient du tenseur de projection par rapport à , le

gradient du tenseur de projection , le gradient du tenseur de
propagation . Or la propagation latérale est réalisée par une addition tensorielle,
on peut donc écrire :

Nous cherchons à calculer , le gradient du tenseur de propagation.

Dans l’annexe B, nous avons vu que le tenseur de propagation, , selon la

direction , et sont sous la formule générale suivante :

Pour la direction :

, est l’ordonnée de l’interface

Pour la direction :

, est l’ordonnée de l’interface

Pour la direction :

, est l’altitude de l’interface

145
Nous pouvons constater que ce tenseur de propagation, , ne dépend pas du point
d’origine du faisceau. Par conséquent, les gradients du tenseur de propagation
par rapport à et sont nuls.

Nous allons maintenant expliciter, pour la direction , le gradient du tenseur de

propagation , désignant les paramètres , et :

Par la symétrie, pour les directions et , il suffit de changer les indices du tenseur et remplacer par
et respectivement.

146
Annexe E Gradient de dose du plan isocentre

Supposons que le plan de traitement possède un isocentre . Les axes de chaque faisceau RCMI se
croisent en ce point. Dans un premier temps, nous considérons que le point isocentre coïncide
avec l’origine du repère fantôme hétérogène. Ceci introduit une contrainte supplémentaire pour le
point d’origine du faisceau RCMI :

Sachant que peut s’exprimer :

Et que vaut :

Nous avons trois contraintes d’égalités :

Soit, après simplification :

(E.1)

La variation de la fonction objectif :

En tenant compte de la contrainte isocentre, les paramètres de l’origine du faisceau ( et ) dépendent

des paramètres d’orientation du faisceau ( et ). Nous avons donc :

147
Par conséquent,

=
Où signifie sous la condition d’isocentre

Maintenant, nous cherchons les quatre dérivées et qui peuvent être

calculées à partir de l’équation (E.1)

La troisième ligne de l’équation (E.1) nous donne :

Alors

Ce qui aux cas particuliers près où (i.e. et valent tout deux 1) donne :

La première ligne de l’équation (E.1) nous donne :

Alors

Nous pouvons en déduire :

148
Annexe F L’évalution d’une distribution de dose par l’index delta

Nous présentons un d’outil de comparaison de doses. La comparaison est réalisée en observant

comment la distribution de dose à évaluer – en général calculée par un algorithme– s’intègre dans une
enveloppe autour de la distribution de référence – obtenue par mesure, par simulation Monte-Carlo ou
comme distribution objectif.

Notons la distribution de référence, et la distribution à évaluer. Une méthode classique de l’état

de l’art est l’index gamma qui compare à en combinant la tolérance en dose (ex. 2mm) et la
tolérance spatiale (ex. 3%). Pour chaque point de référence, une ellipse est formée avec les
dimensions données par et , illustré en Figure 104 , puis pour chaque point de on évalue s’il
est dans l’une des ellipses ou pas.

Figure 104 L’ellipse formée d’un point de référence.

L’index delta, introduit dans (Blanpain, et al., 2009), repose sur le concept de l’enveloppe delta. Cette
enveloppe transforme les tolérances de l’index gamma en dose et en position, en tolérances en dose
uniquement. Cette transformation est faite en calculant, à partir des ellipses de l’index gamma, les
doses minimale et maximale qui peuvent être tolérées en une position donnée (voir Figure 105). Plus
précisément, la dose maximale (resp. minimale) tolérée au point est la dose maximale (resp.
minimale) atteinte par l’ellipse d’un des points de référence situés au voisinage de . L’enveloppe
pour un point donnée est représentée par les flèches dans la Figure 105. Nous noterons (resp.
) la distribution de dose maximale (resp. minimale) tolérée.

Figure 105 L’enveloppe delta.

En chaque point, l’index delta est ensuite calculé à partir de la valeur référence en ce point et des
enveloppes :

149
  Si , alors l’index delta est positif et vaut . Une valeur exacte
de +1 signifie donc que et une valeur exacte de 0 signifie que
.
 Si , alors l’index delta est négatif et vaut . Une valeur exacte
de -1 signifie donc que .

La seconde position sur la Figure 105 est une bonne illustration de l’intérêt de l’index delta vis-à-vis
de l’index gamma classique : dans cette zone de fort gradient, il y a des valeurs de dose qui sont entre
les domaines définies par les ellipses de tolérance. Si la distribution de dose à évaluer a ces valeurs,
une erreur significative sera relevée alors que les valeurs sont logiquement valides. En revanche, ces
valeurs sont comprises dans l’enveloppe delta, ce qui signifierait que ces doses sont belles et bien
valides.

Figure 106 L’enveloppe delta sur un cas de test.

La Figure 106 montre deux distributions de dose 1D : est le build-up sur l’axe d’un faisceau étroit,
et résulte d’un shift de 0.5 mm à droite de puis d’une multiplication par 1.05. L’enveloppe delta
est représentée par deux courbes en pointillé rouge. On remarque immédiatement la praticité du
concept d’enveloppe qui traduit les tolérances et en des tolérances en dose très simples à
visualiser. Des erreurs de calcul (la valeur absolue de l’index delta est supérieure à 1) pour les points
qui sont en dehors de cette enveloppe rouge sont détectées.

Figure 107 Courbes de l’index gamma non interpolé, de l’index gamma interpolé, et de l’index delta.

On voit, sur la Figure 107, les résultats obtenus par trois index : l’index gamma avec et sans
interpolation, ainsi que l’index delta. Ces trois index donnent le même résultat sur la partie droite, où
le gradient est faible. Sur la partie gauche par contre, l’index gamma calculé sans interpolation détecte
des erreurs où il n’y en a pas. L’index delta, aussi bien que l’index gamma interpolé, élimine ces faux
positifs.

150
Les évaluations utilisant une tolérance spatiale (DTA, index gamma, index delta) ont pour objectif
d’autoriser des erreurs un peu plus élevées dans les zones où le gradient est grand. Il a été montré, dans
(Blanpain, et al., 2009), qu’une enveloppe construite à partir de tolérances en dose et en position de
3% et 2 mm respectivement, peut autoriser jusqu’à 25 % d’erreur en dose sur les zones à très fort
gradient, ce qui est très important. Ceci montre toutes les précautions qu’il convient de prendre pour
utiliser des méthodes telles que les index gamma et delta, en particulier lors des évaluations réalisées
en clinique.

151
Bibliographie

AAPM Tissue inhomogeneity corrections for megavoltage photon [Rapport]. - 2004. - 85.
Accuray Cyberknife M6 Series Guide d'implantation [Rapport]. - 2013.
Ahnesjo A Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous
media [Revue] // Med Phys. - 1989. - Vol. 16. - pp. 577-592.
Ahnesjo A et Aspradakis M.M Dose calculations for external photon beams in radiotherapy
[Revue] // Phys Med Biol. - 1999. - Vol. 44. - pp. 99-155.
Ahnesjo A, Andreo P et Brahme A Calculation and application of point spread functions for
treatment planning with high energy photon beams [Revue] // Acta Oncol. - 1987. - Vol. 26. - pp. 49-
56.
Ahnesjo A, Saxner M et Trepp A A pencil beam model for photon dose calculation [Revue] // Med.
Phys. - 1992. - Vol. 19. - pp. 263–273.
ANAES Evaluation de la radiothérapie conformationnelle avec modulation d'intensité [Revue]. -
2003.
Andreo P Monte carlo techniques in medical radiation physics [Revue] // Phys. Med. Biol.. - 1991. -
Vol. 36. - pp. 861–920.
Aventin C [et al.] nouvelles technologies et optimisation des volumes à traiter en radiothérapie
[Revue] // Techniques Hospitalières. - 2012. - pp. 65-72.
Bahi J [et al.] Neural network based algorithm for radiation dose evaluation in heterogeneous
environments [Revue] // Int. Conf. on Artificial Neural Networks. - Athens, Greece : [s.n.], 2006. -
Vol. 4132. - pp. 777–787.
Bey P La radiothérapie sous toutes ses formes. - 2010.
Bjarngard B. E, Rashid H et Obcemea C. H Separation of primary and scatter components of
measured photon beam data [Revue] // Phys. Med. Biol. - 1989. - 12 : Vol. 34 . - pp. 1939-1945.
Bjarngard B.E et Shackford H Attenuation in high-energy x-ray beams [Revue] // Med Phys. - 7 :
Vol. 21. - pp. 1069-1073.
Blake S.W. Artificial neural network modelling of megavoltage photon dose distributions [Revue] //
Phys.Med. Biol. - 2004. - Vol. 49 . - pp. 2515–2526.
Blanpain B et Mercier D The delta envelope: A technique for dose distribution comparison
[Revue] // Med.Phys.. - 2009. - pp. 797-808.
Blanpain B Vers un calcul en temps réel de la dose dans un fantôme segmenté en mailles
homogènes // Thèse : radio physique et imagerie médical. - Toulouse : [s.n.], 2009.
Boyer A et Mok E A photon dose distribution model employing convolution calculations [Revue] //
Med Phys. - 1985. - Vol. 12. - pp. 169-177.
Brahme A Optimized radiation therapy based on radiobiological objectives [Revue] // Semin Radiat
Oncol. - 1999. - 1 : Vol. 9. - pp. 35-47.
Brahme A, Roos JE et Lax I Solution of an integral equation encountered in radiation therapy
[Revue] // Phys.Med.Biol. - 1982. - Vol. 27. - pp. 1221-1229.
Brualla L, Salvat F et Palanco-Zamora R Efficient Monte Carlo simulation of multileaf collimators
using geometry-related varience-reduction techniques [Revue] // Phys Med Biol. - 2009. - 13 : Vol.
54. - pp. 4131-4149.
Burman C [et al.] Planning, delivery, and quality assurance of intensity-modulated radiotherapy
using dunamic multileaf collimator: a strategy for large-scale implementation for the treatment of
carcinoma of the prostate [Revue] // Int J Radiat Oncol Biol Phys. - 1997. - 4 : Vol. 39. - pp. 863-873.
C Martin Benjamin, R Bortfeld Thomas et A Castanon David Accelerating IMRT optimization by
voxel sampling [Revue] // Phys. Med. Biol. - 2007. - pp. 7211-7228.
Childress Nathan [et al.] Dose Calculation Algorithm [Revue] // Mobius3D White Paper. - 2012.
Clarkson J.R A note on depth doses in fields of irregular shape [Revue] // Brit. J. Radiol. - 1941. -
Vol. 14. - pp. 265 - 268.
Colorado Cyberknife Lung Cancer Treatment Colorado Cyberknife [En ligne]. - 2014. -
http://www.coloradocyberknife.com/treatment/lung-and-pulmonary-tumors2/.

152
Convery D.J et Rosenbloom M.E The generation of intensity-modulated fields for conformal
radiotherapy by dynamic collimation [Revue] // Phys Med Biol. - 1992. - 6 : Vol. 37. - pp. 1359-1374.
Coucke P.A Aperçu de la radiothérapie. - 2013.
Craft D Local beam angle optimization with linear programming and gradient search [Revue] // Phys
Med Biol. - 2007. - 7 : Vol. 52. - pp. 127-135.
Cunningham J. R Scatter-air ratios [Revue] // Phys. Med. Biol. - 1972. - 1 : Vol. 17 . - pp. 42-61.
Dieterich S et Gibbs I.C The CyberKnife in clinical use: Current roles, future expectations [Section
du livre] // IMRT,IGRT,SBRT - Advances in the treatment planning and delivery of radiotherapy /
auteur du livre Meyer J.L. - 2011. - Vol. 43.
Dreyfus G [et al.] Réseaux de neurones, Méthodologie et applications [Revue] // Eyrolles. - 2002.
Drzymala R.E [et al.] Dose-Volume Histograms [Revue] // Int.J.Radiation Oncology Biol Phys. -
1991. - 1 : Vol. 21. - pp. 71-78.
Earl M.A [et al.] Inverse planning for intensity-modulated arc therapy using direct aperture
optimization [Revue] // Phys Med Biol. - 2003. - Vol. 48. - pp. 1075-1089.
Emami B et Lyman J Tolerance of normal tissue to therapeutic irradiation [Revue] // Int J Radiat
Oncol Biol Phys. - 1991. - 1 : Vol. 21. - pp. 109-122.
Fogliata A [et al.] On the dosimetric behaviour of photon dose calculation algorithms in the presence
of simple geometric heterogeneities : comparison with monte carlo calculations [Revue] // Phys. Med.
Biol. - 2007. - Vol. 52. - pp. 1363–1385.
Gaede S et Wong E An algorithm for systematic selection of beam directions for IMRT [Revue] //
Med Phys. - 2004. - 2 : Vol. 31. - pp. 376-288.
Gagnon E Le Cyberknife:Radiochirurgie intracrânienne et extracrânienne [Représentation]. - 2012.
Gambini Denus Jean et Granier Robert Manuel pratique de radioprotection [Livre]. - 2007.
Gladwish A [et al.] Segmentation and leaf sequencing for intensity modulated arc therapy [Revue] //
Med Phys. - 2007. - Vol. 34. - pp. 1779-1788.
Habib B [et al.] Evaluation of penfast- a fast Monte carlo code for dose calculations in photon and
electron radiotherapy treatment planning [Revue] // Phys Med. - 2010. - 1 : Vol. 26 . - pp. 17- 25.
Habrand J.-L et Drouet F Tolérance à l’irradiation des tissus sains : moelle épinière [Revue] //
Cancer/Radiothérapie. - 2010. - Vol. 14. - pp. 269–276.
Holder Allen et Salter Bill A tutorial on radiation oncology and optimization [Article] // Mathematics
Faculty Research. - 2004. - 37.
Holder Allen Radiotherapy treatment design and linear programming [Livre]. - [s.l.] : International
series in operations research and management science, 2005. - Vol. 70 : 4 : pp. 741-774.
Horaud R et Monga O Vision par ordinateur [Livre]. - [s.l.] : Editions Hermès, 1995.
ICRU 50 Prescribing, Recording, and Reporting Photon Beam Therapy [Rapport]. - 1993.
ICRU 62 Prescribing, Recording and Reporting Photon Beam Therapy [Rapport]. - 1999.
ICRU 83 Prescribing, Racoding, and Reporting Photon-Beam Intensity-Modulated Radiation Therapy
[Rapport]. - 2010.
IMRT Collaborative Working Group Intensity-Modulated Radiotherapy: Current Status and Issues
of Interest [Revue] // Int J Radiation Oncology Biol Phys. - 2001. - 4 : Vol. 54. - pp. 880-914.
Johns H.E et Cunningham J.R Physics of Radiology [Livre]. - [s.l.] : Thomas books publishing,
1983.
Knoos T [et al.] Limitations of a pencil beam approach to photon dose calculations in lung tissue
[Revue] // Phys. Med. Biol. - 1995. - Vol. 40. - pp. 1411–1420.
Kooy H.M et Rashid H A three-dimensional electron pencil-beam algorithm [Revue] // Phys. Med.
Biol. - 1989. - Vol. 34. - pp. 229–243.
Lafond C Thèse: Analyse et optimisation des performances de la technique VMAT pour son
utilisation en radiothérapie . - Rennes : [s.n.], 2013.
Langen K Tomotherapy's implementation of image-guided adaptive radiaiton therapy. - 2006.
Lefkopoulos D [et al.] La planification inverse en radiothérapie d'intensité modulée [Revue] //
Cancer/Radiother. - 1999. - pp. 160-170.
Löf Johan Development of a general framework ffor optimization of radiation therapy [Livre]. -
2000. - p. 8.
Love P [et al.] The Use of Computers in Radiation Therapy [Livre]. - 2000. - pp. 409-410.

153
Low D.A. [et al.] A technique for the quantitative evaluation of dose distributions [Revue] //
Med.Phys.,25. - 1998. - pp. 656-61.
Mackie Rock Tomotherapy: A new concept for the delivery of dynamic conformal radiotherapy
[Revue] // Med Phys. - 1993. - 6 : Vol. 20.
Mackie T The convolution/superposition method:A model-based dose computation algorithm
[Rapport]. - 2010.
Mackie T, Scrimger J et Battista J A convolution method of calculating dose for 15-mv x rays
[Revue] // Med Phys. - 1985. - Vol. 12. - pp. 188-196.
Mathieu R [et al.] Calculations of dose distributions using a neural network model [Revue] // Phys.
Med. Biol. - 2005. - Vol. 50 . - pp. 1019–1028.
Menguy Y Optimisation quadratique et géométrique de problèmes de dosimétrie inverse // Thèse :
mathématique appliquées. - Grenoble : [s.n.], 1996.
Metcalfe P [et al.] Beam hardening of 10 mv radiotherapy x-rays: analysis using a
convolution/superposition method [Revue] // Phys Med Biol. - 1990. - Vol. 35. - pp. 1533-1549.
Mohan R et Chui C Use of fast fourier transforms in calculating dose distributions for irregularly
shaped fields for three-dimensional treatment planning [Revue] // Med. Phys. - 1987. - Vol. 14. - pp.
70–77.
Mohan R et Wang X The potential and limitations of the inverse radiotherapy technique [Revue] //
Radiother Oncol. - 1994. - 3 : Vol. 32. - pp. 232-248.
Mohan R, Chui C et Lidofsky L Differential pencil beam dose computation model for photons
[Revue] // Med Phys. - 1986. - Vol. 13. - pp. 64-73.
Niemierko A Reporting and analyzing dose distributions: a concept of equivalent uniform dose
[Revue] // Med Phys. - 1997. - 1 : Vol. 24. - pp. 103-110.
O’Connor J. E The density scaling theorem applied to lateral electronic equilibrium [Revue] // Med.
Phys. - 1984. - Vol. 11. - pp. 678–680.
Otto K Volumetric modulated arc therapy: IMRT in a single gantry arc [Revue] // Phys Med. - 2008. -
35. - p. 310.
Ouabri K Développement d'un système de dosimétrie relative des faisceaux de photons de haute
énergie à l'aide de dosimètres thermoluminescents. - 2012.
Pollack A [et al.] Preliminary results of a randomized radiotherapy dose-escalation study comparing
70 Gy with 78 Gy for prostate cancer [Revue] // J Clin Oncol. - 2000. - Vol. 18. - pp. 3904-3911.
Pugachev A [et al.] Role of beam orientation optimization in intensitymodulated radiation therapy
[Revue] // Int J Radiat Oncol Biol Phys. - 2001. - 2 : Vol. 50. - pp. 551-560.
Pugachev A et Xing L Incorporating prior knowledge into beam orientation optimization in IMRT
[Revue] // Int J Radiat Oncol Biol Phys. - 2002. - 5 : Vol. 54. - pp. 1565-1574.
Rowbottom C.G, Nutting C.M et Webb S Beam-orientation optimization of intensitymodulated
radiotherapy: clinical application to parotid gland tumours [Revue] // Radio Oncol. - 2001. - 2 : Vol.
59. - pp. 169-177.
Salvat F [et al.] Penelope-2006, a code system for monte carlo simulation of electron and photon
transport [Revue] // In NEA 6222. - 2006.
Salvat F [et al.] Practical aspects of monte carlo simulation of charged particle transport : Mixed
algorithms and variance reduction techniques [Revue] // Radiat. Environ. Biophys.. - 1998. - Vol. 38. -
pp. 15-22.
Schweikard A, Schlaefer A et Adler J.R Resampling: A optimization method for inverse planning in
robotic radiosurgery [Revue] // Med Phys. - 2006. - Vol. 33. - p. 4005.
SFPM Contrôles de qualité en radiothérapie conformationnelle avec modulation d'intensité
[Rapport]. - 2009.
Shepard D.M [et al.] An arc-sequencing algorithm for intensity modulated arc therapy [Revue] //
Med Phys. - 2007. - Vol. 34. - pp. 464-470.
Shepard David M [et al.] Optimizing the delivery of radiation therapy to cancer patients [Revue] //
Society for Industrial and Applied Mathematics Review. - 1999. - 4 : Vol. 41. - pp. 721-744.
Soderstrom S et Brahme A Which is the most suitable number of photon beam portals in coplanar
radiation therapy? [Revue] // Int J Radiat Oncol Biol Phys. - 1995. - 1 : Vol. 33. - pp. 151-159.
Spirou SV et Chui CS A gradient inverse planning algorithm with dose-volume constraints [Revue] //
Med Phys. - 1998. - 3 : Vol. 25. - pp. 321-333.

154
Stein J [et al.] Number and orientations of beams in intensity-modulated radiation treatments
[Revue] // Med Phys. - 1997. - 2 : Vol. 24. - pp. 149-160.
Tillikainen L [et al.] A 3d pencil-beam-based superposition algorithm for photon dose calculation in
heterogeneous media [Revue] // Phys. Med. Biol. - 2008. - Vol. 53. - pp. 3821–3839.
Tombropoulos R.Z, Adler J.R et Latombe J.C CARABEAMER: a treatment planner for a robotic
radiosurgical system with general kinematics [Revue] // Medical Image Analysis. - 1999. - 3 : Vol. 3. -
pp. 237-264.
TomoTherapy Tomotherapy - we believe in a better way [En ligne]. - 2011. - 28 janvier 2014. -
www.tomotherapy.com/video/tags/tag/tomotherapy.
Ulmer W, Pyyry J et Kaissl W A 3d photon superposition/ convolution algorithm and its foundation
on results of monte carlo calculations [Revue] // Phys. Med. Biol. - 2005. - Vol. 50. - pp. 1767–1790.
Van Dyk J [et al.] Tomotherapy: A " revolution" in radiation therapy [Revue] // Physics in Canada. -
2002.
Vanderstraeten B, Reynaert N et Paelinck L Accuracy of patient dose calculation for lung imrt : A
comparison of monte carlo, convolution/superposition, and pencil beam computations [Revue] // Med.
Phys. - 2006. - Vol. 33. - pp. 3149–3158.
Venselaar J, Welleweerd H et Mijnheer B Tolerances for the accuracy of photon beam dose
calculations of treatment planning systems [Revue] // Radio and Oncol. - 2001. - Vol. 60. - pp. 191-
201.
Vuillez Jean-Philippe Biophysique [Livre]. - 2009. - Vol. chapitre 4 Dosimétrie.
Wang XH et Mohan R Optimization of intensity modulated 3D conformal treatment plans based on
biological indices [Revue] // Radiother Oncol. - 1995. - 2 : Vol. 37. - pp. 140-152.
Wong E, Chen J et Greenland J Intensity-modulated arc therapy simplified [Revue] // Int J Radiat
Oncol Biol Phys. - 2002. - Vol. 53. - pp. 222-235.
Wu X et Zhu Y A neural network regression model for relative dose computation [Revue] // Phys.
Med.Biol.. - 2000. - Vol. 45 . - pp. 913–922 .
Yu C.X [et al.] Clinical implementation of intensity-modulated arc therapy [Revue] // Int J Radiat
Oncol Biol Phys. - 2002. - Vol. 53. - pp. 453-463.
Yu C.X Intensity-modulated arc therapy with dynamic multileaf collimation: An alternative to
tomatherapy [Revue] // Phys Med Biol. - 1995. - Vol. 40. - pp. 1435-1449.

155
Publications et communications
Brevet
Yang M., Mercier D., Bordy J.M.: Procédé d'optimisation d'une planification de traitement en
radiothérapie. Brevet CEA.

Conférence
Yang M., Blanpain B., Mercier D., Bordy J.M.: A new method for real time dose calculation with
partition of the 3D phantom into homogeneous meshes. Présentation poster à European Society for
Radiotherapy and Oncology (ESTRO), Barcelone, Espagne, 9-13 Mai, 2012.

Séminaire
Présentation orale et posters aux journées des doctorants de l’école doctorale MIPEGE en mars 2010,
2011 et 2012.

156