L Pibouleau
email: luc.pibouleau@ensiacet.fr

S Domenech

A Davin

C Azzaro-Pantel

Expérimentations numériques sur les variantes et paramètres de la méthode du recuit simulé

Keywords: Recuit simulé, Implémentation, Codage, Paramètres de contrôle, Plan d'expériences Simulated annealing, Implementation, Code, Control parameters, Design of experiments

 sont les variantes les plus efficaces. Parmi les trois codages testés, aucun ne supplante les autres. La température initiale, la longueur du palier de température, le facteur multiplicatif de décroissance de la température et le test d'arrêt ont été étudiés au moyen d'un plan factoriel 2 4 . Le paramètre le plus important est la longueur du palier de température, le moins important est la température initiale, et le coefficient de décroissance de la température et le nombre de paliers pour effectuer l'arrêt ont une importance similaire.

Introduction

Contrairement aux méthodes classiques de la recherche opérationnelle (programmation mathématique, programmation mixte, programmation dynamique, méthodes arbores-centes), les algorithmes stochastiques de type Recuit Simulé (RS) n'exigent aucune connaissance des propriétés mathématiques du problème; il est seulement nécessaire de pouvoir évaluer la qualité d'une solution. Il s'agit là d'un point qui plaide en leur faveur pour des applications en Génie des Procédés, où les fonctions objectifs sont généralement évaluées à partir de modèles de simulation complexes, dont aucune expression analytique n'est Nomenclature A coefficient utilisé pour calculer la valeur initiale de la température B coefficient utilisé pour calculer la longueur du palier de température C coefficient utilisé dans le schéma de décroissance de la température D nombre de paliers successifs sur lesquels la solution courante demeure inchangée (test d'arrêt) ddl nombre de degrés de liberté f fonction objectif f * valeur de la fonction objectif pour la solution théorique fa valeur minimale de la fonction objectif déterminée par la procédure du recuit simulé GA probabilité de [START_REF] Gauber | Time-dependent statistics of the Ising Model[END_REF], décroissance de la température de [1] GK probabilité de [START_REF] Gauber | Time-dependent statistics of the Ising Model[END_REF], décroissance de la température de [START_REF] Kirkpatrick | Simulated annealing: theory and applications[END_REF] h amplitude de l'intervalle de confiance pour sm2 H matrice d'Hadamard I matrice identité MA probabilité de [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], décroissance de la température de [1] MK probabilité de [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], décroissance de la température de [ disponible, ce qui rend impossible l'étude de leurs propriétés mathématiques. De plus, ces algorithmes sont également faciles à mettre en oeuvre (le moteur est commun, il y a moins de programmation spécifique à chaque problème traité). Le Tableau 1 récapitule les différences entre une approche de type recherche opérationnelle et la résolution du même problème par un algorithme stochastique. Ce Tableau semble plaider de fac ¸on irrévocable en faveur des algorithmes stochastiques, mais le fait que la précision des résultats ne soit pas maîtrisée, car les conditions d'arrêt ne le sont pas, constitue un handicap pour ces derniers. Par ailleurs, les algorithmes stochastiques comportent un grand nombre de paramètres, dont un mauvais réglage peut conduire à des solutions de piètre qualité. Néanmoins, ils sont très efficaces pour la détermination rapide de bonnes solutions à un problème posé, c'est pourquoi ils sont actuellement très en vogue.

Applications en génie des procédés

Au cours des deux dernières décennies, les méthodes d'optimisation stochastiques ont largement été mises en oeuvre dans divers domaines du génie des procédés, comme le montre l'étude statistique effectuée par [START_REF] Pibouleau | A survey of optimisation tools through ESCAPE symposia[END_REF]. Sans vouloir être exhaustif, il existe à titre d'exemple, des applications à la conception de séquences de séparation [START_REF] Floquet | Separation sequence synthesis: how to use simulated annealing procedure?[END_REF][START_REF] Laquerbe | Conception de séquences de séparation par un algorithme génétique[END_REF], de réseaux d'échangeurs de chaleur [START_REF] Dolan | Algorithmic efficiency of simulated annealing for heat exchanger network design[END_REF][START_REF] Athier | Synthesis of optimum heat exhanger networks by simulated annealing[END_REF], de réacteurs pour l'industrie chimique [START_REF] Cordero | Conception optimale de réseaux de réacteurs par une méthode mixte[END_REF], de distribution d'énergie [START_REF] Dolan | Process optimization via simulated annealing: application to network design[END_REF], aux problèmes d'ordonnancement [START_REF] Peyrol | Simulation par événements discrets d'un atelier de fabrication de composants électroniques[END_REF], ou à la modélisation des écoulements dans des locaux ventilés [START_REF] Laquerbe | Synthesis of RTD Models via stochastic procedures: simulated annealing and genetic algorithms[END_REF]. D'autres applications concernent l'optimisation [START_REF] Baudet | Ordonnancement à court terme d'un atelier discontinu de chimie fine[END_REF], la conception [START_REF] Bernal-Haro | Multi-objective batch plant design: a two-stage methodology. Part II. Development of a genetic algorithm and results analysis[END_REF] ou le remodelage [START_REF] Dedieu | A retrofit design strategy for multipurpose batch plants[END_REF] d'ateliers discontinus de chimie fine. Le codage des solutions potentielles d'un problème donné constitue un point d'autant plus important des méthodes d'optimisation stochastiques, que le nombre de contraintes est élevé, ce qui est souvent le cas dans les applications en génie des procédés. La stratégie utilisée dans les exemples précédents consistait à introduire directement les contraintes dans le codage. Toutes ces études font ressortir la difficulté de choisir, d'une part entre les variantes possibles, et d'autre part de définir des valeurs des paramètres de contrôle pour obtenir des solutions de bonne qualité [START_REF] Davin | Optimisation des paramètres du recuit simulé en génie des procédés appliquée aux réseaux d'appareillages (réacteurs, échangeurs de chaleur)[END_REF].

Méthode du recuit simulé

Principes de la méthode

La méthode du recuit simulé, popularisée par [START_REF] Kirkpatrick | Simulated annealing: theory and applications[END_REF], est une procédure d'optimisation multi-variables développée depuis une vingtaine d'années environ. Contrairement aux méthodes déterministes, qui sont pour la plupart des méthodes locales, il s'agit d'une procédure évolutive qui permet d'éviter de se bloquer dans des zones de minima locaux, en autorisant occasionnellement un accroissement de la fonction objectif. Elle se fonde sur une analogie entre les problèmes d'optimisation à caractère combinatoire et ceux de la physique statistique représentant le comportement d'une opération thermodynamique, le recuit physique des solides.

Le recuit physique des métaux

Considérons par exemple le cas de la croissance d'un monocristal. La technique du recuit consiste à chauffer préalablement le matériau pour le porter quasiment à l'état liquide. L'état thermodynamique du système peut être traduit par une fonction thermodynamique: l'énergie libre. A l'état liquide, cette énergie est élevée.

Le refroidissement lent du système, en marquant des paliers de température, redistribue l'arrangement des atomes de la matière jusqu'à atteindre un état énergétique faible, celui de grands cristaux, soit un état stable. Cet état stable correspond au minimum absolu de l'énergie libre (E). Si la vitesse de refroidissement n'est pas contrôlée correctement, ou si la température initiale n'est pas suffisamment élevée, le solide n'atteint pas l'équilibre thermique à chaque palier de température; dans de telles circonstances, il est obtenu une structure amorphe, état métastable qui correspond à un minimum local de l'énergie libre. Il est donc nécessaire d'atteindre l'équilibre thermodynamique pour chaque palier de température, qui est caractérisé par une fonction de distribution de Boltzmann des états énergétiques. La probabilité pour un système de posséder une énergie donnée E à la température T est ainsi proportionnelle à exp(-E/k b T ) où k b est la constante de Boltzmann. A proximité du zéro absolu, seuls les états de faible énergie ont une probabilité non nulle d'exister.

Ainsi même à basse température, une transition peut subvenir d'un niveau énergétique bas vers un niveau plus haut. La probabilité de telles transitions est d'autant plus faible que le niveau de température est bas. Ce sont ces transitions qui sont supposées être à l'origine des mécanismes permettant d'atteindre in fine un état d'énergie minimale au lieu de conduire à un état métastable.

Transposition du recuit simulé comme technique d'optimisation

A partir d'un arrangement moléculaire initial (dans ce cas, une solution acceptable S du problème d'optimisation), d'état énergétique E i donné correspondant à une température T i particulière, un nouvel arrangement moléculaire (solution acceptable S ) possédant une énergie E j est généré de manière stochastique. Si la nouvelle solution acceptable a pour effet d'abaisser l'énergie, E j ≤ E i (c'est à dire de diminuer le critère, f(S ) < f(S), le changement s'opère car l'état est plus stable. Si, au contraire, l'état énergétique est plus élevé, E j ≥ E i (le critère croît), la transition peut être acceptée ou non selon une règle statistique. Ces règles d'acceptation permettent de s'extraire de zones à minima locaux, en effet:

• à température élevée, cette probabilité est grande, c'est-àdire que de très nombreuses solutions dégradant le critère seront conservées, on assure ainsi une large exploration de l'espace de recherche; • à température basse, cette probabilité tend vers zéro, donc pratiquement toutes les solutions dégradant le critère seront refusés; la procédure tend alors à se comporter comme une méthode de descente classique.

L'implémentation de la procédure du recuit simulé nécessite de fixer quatre paramètres, qui jouent un rôle déterminant sur son efficacité.

(i) La température initiale. (ii) La longueur du palier de température. (iii) Le schéma de décroissance de la température. (iv) Le nombre de paliers successifs sur lesquels la solution courante S n'évolue pas, et provoque l'arrêt de la recherche.

Objectifs de l'étude

A partir de fonctions mathématiques dont les solutions sont connues, des essais numériques systématiques, fondés sur les plans d'expériences, ont été effectués afin de dégager des idées sur la fac ¸on de fixer les divers paramètres, ainsi que les variantes de mise en oeuvre de la méthode du recuit simulé. En essayant de dégager des grandes orientations dans le choix des paramètres de cet algorithme, cette étude vise à pallier l'empirisme qui est généralement présent dans la mise en oeuvre de ces méthodes, et qui peut conduire à un très grand nombre d'essais avant de déterminer un jeu de paramètres adéquats.

Variantes de mise en oeuvre et paramètres de la procédure du recuit simulé

A partir d'une solution acceptable S, codée en binaire, du problème d'optimisation, une nouvelle solution acceptable S est généré de manière stochastique. Si la nouvelle solution acceptable a pour effet de diminuer le critère, f(S ) < f(S), le changement s'opère car l'état est plus stable. Si, au contraire le critère croît, la nouvelle solution S peut être acceptée ou non selon une règle statistique. La bibliographie fait ressortir l'existence de deux protocoles classiques communément utilisés:

• la procédure de [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], notée , où l'acceptation de la nouvelle solution est réalisée avec une probabilité égale à:

( f = (f (S ) -f (S))) = exp -f T i si f ≥ 0; p = 1 si f < 0 (1)
• la variante proposée par [START_REF] Gauber | Time-dependent statistics of the Ising Model[END_REF], notée , où la probabilité d'acceptation est uniforme, et définie par:

p( f = (f (S ) -f (S))) = exp -f T i /(1 + exp(-f/T i )) (2)
En appliquant une de ces règle d'acceptation, qui permettent de s'extraire de zones à minima locaux, sur un palier de température constant, il est assuré que l'évolution du système vers son équilibre thermodynamique, correspond au minimum d'énergie libre atteignable pour la valeur de la température considérée. Lorsque cet équilibre est atteint, après un nombre "suffisant" de transformations élémentaires, la température est abaissée "légèrement" selon un schéma de décroissance donné, le systéme évolue à nouveau pour aboutir à un nouvel état d'équilibre. Il existe deux schémas classiques de décroissance de la température:

• la décroissance géométrique, proposée par [START_REF] Kirkpatrick | Simulated annealing: theory and applications[END_REF], notée :

T i+1 = CT i avec 0 < C < 1 (3)
• la procédure utilisée par [START_REF] Aarts | Statistical cooling: a general approach to combinatorial optimization problems[END_REF], notée :

T i+1 = T i (1 + (T i ln(1 + C)/3σ)) (4)
où σ désigne l'écart-type de l'ensemble des valeurs calculées de la fonction objectif à la température courante et C, le paramètre qui contrôle la vitesse du recuit; plus C est faible, plus le recuit est lent. Ces procédures sont discutées par [START_REF] Dolan | Algorithmic efficiency of simulated annealing for heat exchanger network design[END_REF][START_REF] Patel | Preliminary design of multiproduct noncontinuous plants using simulated annealing[END_REF].

Outre les variantes qui viennent d'être décrites, l'implémentation de la procédure du recuit simulé nécessite de fixer quatre paramètres, qui jouent un rôle déterminant sur son efficacité.

• La température initiale, dont la dimension (en unité) est celle de la fonction objectif, et dont la valeur dépend directement de la valeur initiale de cette fonction. Son choix devra être tel que la plupart des dégradations sont autorisées en début de procédure, et que seules les transitions améliorant le critère sont acceptées en fin de procédure. Une étude [START_REF] Athier | Synthesis of optimum heat exhanger networks by simulated annealing[END_REF] a montré qu'une valeur initiale de la température de l'ordre de celle de la fonction objectif initiale, conduit à de bons résultats. Ainsi, l'initialisation proposée est définie par:

T = f (S) A ( 5
)
où f(S) représente la valeur du critère pour la solution initiale S, et A un paramètre à fixer. • La longueur du palier, qui représente le nombre de solutions testées avant d'appliquer le schéma de décroissance de la température, est supposée constante pour l'ensemble des niveaux de température. Si cette longueur est trop faible, l'exploration de l'espace de recherche peut être trop partielle, alors qu'une valeur trop importante peut avoir pour effet de ralentir, voire même de bloquer la recherche. La stratégie proposée consiste à choisir la longueur du palier proportionnelle à la longueur du codage de l'ensemble des n variables du problème:

L = B n i=1
(Longueur du codage de la variable i) (6)

B est le deuxième paramètre à choisir.

• Le troisième paramètre est le coefficient C de réduction de la température, qui doit être tel que:

0 < C < 1 (7)
Une valeur trop faible pour ce coefficient peut provoquer une exploration très partielle de l'espace de recherche. Les valeurs classiques de ce coefficient se situent dans l'intervalle [0.75, 0.95]. • La procédure du recuit simulé est souvent arrêtée lorsque la température devient inférieure à une valeur donnée. Le choix de cette valeur limite, qui conditionne l'efficacité de la recherche, est très délicat car ce terme est étroitement lié aux autres paramètres de la méthode (température initiale, longueur des paliers, schéma de réduction de la température). Ainsi, un test d'arrêt différent portant sur la non évolution de la solution courante S sur D paliers consécutifs a été mis en oeuvre; ce quatrième paramètre D a été choisi tel que:

0 < D = 9 (8)
L'algorithme de la méthode du recuit simulé est défini par les étapes suivantes.

L'objectif de ce travail consiste en premier lieu à étudier sur des fonctions mathématiques dont les solutions sont connues (fonctions utilisées par [START_REF] Alliot | Intelligence Artificielle et Informatique Théorique[END_REF]), les diverses variantes de la méthode du recuit simulé, reportées ci-dessous:

• GA: Probabilité de [START_REF] Gauber | Time-dependent statistics of the Ising Model[END_REF], décroissance de la température de [START_REF] Aarts | Statistical cooling: a general approach to combinatorial optimization problems[END_REF] • GK: Probabilité de [START_REF] Gauber | Time-dependent statistics of the Ising Model[END_REF], décroissance de la température de [START_REF] Kirkpatrick | Simulated annealing: theory and applications[END_REF] • MA: Probabilité de [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], décroissance de la température de [START_REF] Aarts | Statistical cooling: a general approach to combinatorial optimization problems[END_REF] • MK: Probabilité de [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], décroissance de la température de [START_REF] Kirkpatrick | Simulated annealing: theory and applications[END_REF] Ensuite les effets et interactions des quatre paramètres de contrôle (A, B, C et D) seront analysés en utilisant la technique des plans d'expériences centrés.

Fonctions étudiées et divers codages

Pour chaque cas traité, la longueur du codage dépend conjointement de l'intervalle de variation des variables, c'est à dire de l'espace de recherche, et de la précision requise sur ces variables.

Dans toute la suite, il sera considéré que la convergence est atteinte lorsque |f * -fa| < 10 -3 où f * et fa désignent respectivement la solution théorique et la solution approchée par la procédure du recuit simulé.

Les quatre premières fonctions étudiées possèdent plusieurs minima locaux, alors que les deux dernières comportent uniquement un optimum global, mais sont délicates à optimiser. Pour chacune de ces fonctions, trois codages arbitraires sont utilisés, afin d'analyser leur influence. Le codage 1, dans lequel le signe et la partie entière sont codés en binaire, et la première décimale est représentée par quatre bits de poids respectifs 0.1; 0.2; 0.4 et 0.6, dont la somme est supérieure à l'unité. Pour la deuxième décimale, les poids sont divisés par 10, et ainsi de suite. Le codage 2 est analogue au codage 1 en ce qui concerne le signe et la partie entière, et la première décimale est codée par quatre bits de poids respectifs 0.1, 0.2, 0.3 et 0.3, dont la somme vaut 0.9, mais dans laquelle le terme 0.3 apparaît deux fois, ce qui peut créer un biais. Pour la deuxième décimale, les poids sont divisés par 10, et ainsi de suite. Enfin, le codage 3, qui est également analogue au codage 1 pour son signe et sa partie entière, représente la partie décimale par un codage binaire de type entier, qui est par exemple divisé par 1000, si l'on désire coder trois décimales. Les valeurs réelles des variables sont utilisées pour calculer la fonction objectif. La procédure de calcul des valeurs réelles des variables est analogue pour tous les autres codes et toutes les autres fonctions.

Fonction dite "boîte d'oeufs"

Nom: BOITE f (x1, x2) = -x1 sin( |x1|) -x2 sin( |x2|) x1 ∈ [-512, 512], x2

Fonction de Griewank

Nom: GRIE

f (x 1 , x 2 ) = 0.00025(x1 2 + x2 2 ) -cos(x1) cos x2 √ 2 x1 ∈ [-π, π], x2 ∈ [-π, π] f * = -1 x1 * = x2 = 0
où les composantes x1 et x2 sont définies à 10 -7 près.

Pour les codages 1 et 2, la longueur du codage pour une variable est égale à 34; pour le codage 3, cette longueur vaut 31 (Fig. 2). Pour les codages 1 et 2, la longueur du codage pour une variable est égale à 20; pour le codage 3, cette longueur vaut 18 (Fig. 3).

Fonction de De Jong

Fonction de Michalevitch

Nom: MICHA

f (x) = -sin(x) sin 10x 2 π 20 x ∈ [0, π] f * = -1 x * = π 2
où la composante x est définie à 10 -7 près. La variable étant positive, le signe n'est pas codé; pour les codages 1 et 2, la longueur du codage est égale à 33; pour le codage 3, cette longueur vaut 30 (Fig. 4).

Fonction de Pibouleau

Nom: PIB

f (x1, x2) = (1 -x1) 2 + k(x2 2 -x1 2 ) aveck = 0.001 x1 ∈ [-π, π], x2 ∈ [-π, π] f * = - k 1 -k = -0.001 x1 * = 1 1 -k , x2 * = 0
où les composantes x1 et x2 sont définies à 10 -7 près.

Pour les codages 1 et 2, la longueur du codage pour une variable est égale à 34; pour le codage 3, cette longueur vaut 31 (Fig. 5).

Fonction de Rosenbrock

Nom: ROS

f (x1, x2) = 100(x1 2 -x2) 2 + (1 -x1) 2 x1 ∈ [-π, π], x2 ∈ [-π, π] f * = 0 x1 * = x2 * = 1
où les composantes x1 et x2 sont définies à 10 -7 près.

Pour les codages 1 et 2, la longueur du codage pour une variable est égale à 34; pour le codage 3, cette longueur vaut 31 (Fig. 6).

Génération d'une solution S voisine de la solution courante S

Les problèmes testés étant purement mathématiques, il est supposé que les codages de chaque variable sont identiques pour toutes, ce qui peut ne pas être le cas pour des problèmes physiques où les variables ont des ordres de grandeurs différents. Ainsi le vecteur binaire S représentant l'ensemble des variables du problème est de longueur m avec m = n × lcode où n et lcode désignent respectivement le nombre de variables du problème et la longueur du codage de chaque variable.

La génération d'une solution S voisine de la solution courante S est effectuée comme suit. Un nombre entier k compris entre 1 et m est généré selon une loi uniforme; pour définir S , on remplace la k iéme composante du vecteur S par son complément binaire.

Comparaison des méthodes

Essais effectués

Pour chaque fonction présentée ci-dessus, les valeurs choisies pour les quatre paramètres A, B, C et D sont reportées dans le Tableau 2.

En utilisant la technique des plans d'expériences centrés, il existe 17 combinaisons de paramètres (en tenant compte du centre des essais présenté en dernière position), comme il est mentionné dans le tableau ci-dessous. Pour chaque combinaison de ces paramètres, cinq essais consécutifs avec initialisations aléatoires ont été réalisés, ce qui conduit à 85 essais par fonction (Tableau 3).

Etude des succès

Rappelons que par définition, un succès est tel que |f * -fa| < 10 -3 où f * et fa désignent respectivement la solution théorique et la solution approchée par la procédure du recuit simulé.

Dans les tableaux de résultats, les conventions suivantes ont été utilisées (Tableaus 4-6):

• Valeurs en caractères normaux: pourcentage de succès (resp. échecs) sur 85 essais pour chaque fonction. • Moy-Méth: pourcentage moyen de succès (resp. échecs) sur les six fonctions, pour chaque méthode. • Moy-Fct: pourcentage moyen de succès (resp. échecs) sur les quatre méthodes, pour chaque fonction.

Etude des échecs notoires

Par définition, un échec notoire est tel que |f * -fa| > 1 où f * et fa désignent respectivement la solution théorique et la solution approchée par la procédure du recuit simulé (Tableaus 7-9).

Discussion

De l'étude des succès, il peut être conclu que, quelque soit le codage:

• la méthode MK est la plus performante;

• la méthode de [START_REF] Aarts | Statistical cooling: a general approach to combinatorial optimization problems[END_REF] utilisée pour faire évoluer la température est la moins efficace.

Ces deux conclusions sont confirmées par l'étude des échecs.

Il est évident que pour traiter un problème physique, le codage est un élément primordial, qui doit pouvoir représenter sans ambiguïté et de manière biunivoque l'ensemble des solutions acceptables. Dans le cas présent, les trois codages mis en oeuvre satisfont ce principe, ont des longueurs comparables, et l'on ne note pas de différence significative sur les pourcentages moyens de succès et d'échecs par fonction selon le codage utilisé.

Par la suite l'étude sera poursuivie avec la méthode MK, en utilisant par exemple le codage 1.

Etude des effets et des interactions

Données de l'étude

Pour étudier les effets des quatre paramètres cités cidessus, ainsi que leurs interactions, un plan factoriel complet à quatre facteurs a été mis en oeuvre; son format général est donné dans le Tableau 10, où la sous-matrice carrée d'ordre 16 définie par les effets, les interactions et la moyenne constitue une matrice d'Hadamard, notée H, qui vérifie donc:

H T H = 16I (9)
Pour chaque combinaison de paramètres, le terme sm1 (second membre 1) représente le pourcentage de succès sur les cinq essais réalisés avec initialisations aléatoires, et le terme sm2 donne le nombre moyen (exprimé en milliers) d'évaluations de la fonction objectif. La valeur de sm2 n'est pas reportée lorsque sm1 est nul. Les valeurs des seconds membres sont reportées dans le tableau ci-après pour chacune des six fonctions étudiées (Tableau 11).

Effets et interactions

Les effets et interactions, pour le critère sm1, puis pour le critère sm1/sm2, des variables du plan d'expériences sont reportées dans les tableaux ci-dessous (Tableau 12).

En considérant le seul critère sm1 (pourcentage de succès), on ne peut pas dégager de conclusion nette sur les effets moyens des paramètres A, B, C et D; il en est de même pour leurs interactions (Tableau 12). Ce constat conduit à utiliser un second critère, constitué par le rapport sm1/sm2 (terme supposé nul lorsque sm1 est nul), qui semble mieux représenter l'efficacité d'une méthode (Tableau 13).

Analyse qualitative

Le classement des effets et interactions pour le critère sm1/sm2 est reporté dans le tableau ci-dessous (Tableau 14

B-D-C-A CD BD AB AC AD BC BCD ABD ABC ACD De Jong B-C-D-A BC BD CD AC AD AB BCD ABD ABC ACD Grie B-C-D-A BC BD CD AB AC AD BCD ABC ABD ACD Micha B-D-C-A CD AD BD BC AB AC ABC ACD ABD BCD Pib B-C-D-A BC BD CD AD AC AB BCD ABD ABC ACD Ros B-C-D-A AD BC CD BD AB AC ACD ABD BCD ABC

il semble normal que le paramètre A ait peu d'influence sur la convergence. Ce constat est confirmée par les faibles interactions d'ordre deux faisant intervenir A. Après l'élimination de ce paramètre, on constate, d'après les effets moyens et les interactions, que les paramètres C et D se comportent sensiblement de la même manière, l'effet moyen du paramètre C étant toutefois légèrement plus important.

Analyse quantitative: essais au centre du plan (A = 1, B = 1, C = 0.85, D = 5)

Dans la mesure où pour chaque fonction, au plus cinq mesures par point du plan sont disponibles pour le terme sm1/sm2 (en cas de non convergence pour un essai, le nombre d'évaluations de la fonction objectif n'est pas connu), le test de Fisher portant sur l'égalité des variances est non significatif, et il peut donc être considéré que la variance des mesures σ2(sm2), demeure constante pour l'ensemble des points du plan. La variance commune σ 2 de sm2 pour les effets moyens, les interactions, et la moyenne calculées à partir des 16 mesures, est alors donnée par: Le nombre de réponses assurant la convergence variant selon l'exemple, l'analyse est effectuée fonction par fonction. Dans un premier temps cinq essais (c'est à dire cinq exécutions de la procédure du recuit simulé, avec des initialisations aléatoires) par fonction ont été effectués, et pour certaines (Boîte, Micha, Pib), le nombre d'essais concluants était égal à un (conduisant à un nombre de degrés de liberté (ddl) infini), ou à deux (nombre de ddl égal à un, fournissant une valeur très élevée dans la table de Student). Pour pallier cette difficulté, le nombre d'essais a été augmenté à 10.

σ 2 = 1 16 (16σ 2 (sm2)) = 1 16 (σ 2 (sm2)) (10a) sachant que: sm2 = 1 n n i=1 nev i ; ddl = n -1 (10b) σ 2 (sm2) = 1 n(n -1) n i=1 (nev i -sm2) 2 (
La variable considérée est le rapport sm1/sm2, où l'intervalle de confiance pour sm2 a une amplitude égale à h; l'intervalle de confiance pour le rapport sm1/sm2 est donné par:

sm1 sm2 + h , sm1 sm2 -h ( 11
)
Cet intervalle de confiance n'est pas symétrique par rapport à la valeur de sm1/sm2. On peut aisément se ramener au cas classique d'un intervalle symétrique, en supposant l'écart-type de sm2 faible devant sm2, hypothèse qui a été confirmée par les divers calculs. Soit α1 le terme:

α1 = sm1 sm2 + h = sm1 1 sm2 × 1 1 + h/sm2 ≈ sm1 1 sm2 × 1 - h sm2 (12)
il vient:

α1 = sm1 1 sm2 - h sm 2 2 = sm1 sm2 - h × sm1 sm 2 2 = sm1 sm2 -β (13)
avec:

β = h × sm1 sm2 2 (14)
De même, on a:

α2 = sm1 sm2 -h = sm1 sm2 + β ( 15
)
où β est toujours défini par la relation [START_REF] Laquerbe | Conception de séquences de séparation par un algorithme génétique[END_REF]. Les valeurs de sm1 et sm2 sont mentionnées dans le Tableau 15. Les nombres d'évaluations (exprimés en milliers)

Conclusion

Les principales variantes de mise en oeuvre de la méthode du recuit simulé, divers codages des variables, ainsi que les quatre principaux paramètres de contrôle (température initiale, longueur du palier de température, décroissance de la température et test d'arrêt) ont été étudiés à partir de six fonctions dont les solutions analytiques sont connues.

Un plan factoriel complet 2 4 a été mis en oeuvre pour définir les variations des paramètres de contrôle, et pour Dans un premier temps, l'étude du taux de convergence a permis de constater qu'au niveau de la mise en oeuvre, la procédure de [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] couplée avec le schéma de décroissance de la température de [START_REF] Kirkpatrick | Simulated annealing: theory and applications[END_REF] fournissait les meilleurs résultats, alors que le schéma de décroissance de la température de [START_REF] Aarts | Statistical cooling: a general approach to combinatorial optimization problems[END_REF] conduisait aux plus mauvais. Cette conclusion a été confortée par l'étude du taux d'échecs (c'est à dire de non convergence).

Trois codages pour les variables ont été également utilisés, mais aucune conclusion nette au sujet de leur influence sur le taux de convergence n'a pu être dégagée. Il est évident que pour traiter un problème physique, le codage est un élément primordial, qui doit pouvoir représenter sans ambiguïté et de manière biunivoque l'ensemble des solutions acceptables. Dans le cas présent, les trois codages mis en oeuvre satisfont ce principe, ont des longueurs comparables, et aucune différence significative n'est notée sur les pourcentages moyens de succès et d'échecs par fonction selon le codage utilisé.

Ces deux constats ont conduit à poursuivre l'étude des effets et interactions des quatre paramètres de contrôle avec un seul codage et l'implémentation de Metropolis -Kirkpatrick. Le plan factoriel complet 2 4 a été mis en oeuvre pour les six fonctions avec cinq essais à chaque sommet du plan, et avec deux types de réponses: le pourcentage de succès (sm1) et l'efficacité définie par le rapport sm1/sm2 où sm2 représente le nombre d'évaluations (exprimé en milliers) de la fonction objectif. En considérant la seule réponse sm1, aucune conclusion nette n'a pu être dégagée sur les effets moyens et interactions des paramètres, alors que la réponse sm1/sm2 a permis de constater qualitativement et quantitativement des tendances claires. L'étude quantitative a été réalisée pour chaque fonction, de fac ¸on classique au centre du plan à partir de 10 essais avec initialisations aléatoires.

Le calcul des intervalles de confiance a permis de déterminer les seuils de signification des effets et interactions. Pour l'ensemble des six fonctions, l'effet moyen le plus important est dû au paramètre B (qui définit la longueur du palier de température en fonction de la longueur du codage des variables), tandis que le paramètre A (qui permet de calculer la température initiale en fonction de la valeur initiale de la fonction objectif) avait l'effet moyen le plus faible. En effet, le paramètre A donne la température initiale, et dans la mesure où le test d'arrêt est effectué sur la non évolution de la fonction objectif sur D paliers consécutifs, et non sur une valeur minimale de la température, il semble normal que ce paramètre ait peu d'influence sur la convergence. Ce constat a été confirmée par la non signification des interactions d'ordre deux faisant intervenir A, ainsi que par la non signification de l'interaction d'ordre quatre. Les paramètres C (coefficient de décroissance géométrique de la température) et D se comportent sensiblement de la même manière, avec un léger avantage pour C.

En conclusion, et en l'absence de toute autre information, il peut être préconisé d'utiliser la méthode du recuit simulé avec les éléments suivants:

• acceptation d'une solution dégradant la fonction objectif: règle de Metropolis; • décroissance de la température: schéma géométrique de Kirkpatrick; • valeur initiale de la température: valeur de la fonction objectif pour la solution initiale (A = 1); • longueur du palier de temperature: somme des longueurs des codages de chaque variable (B = 1); • coefficient de réduction de la température: 0.85 (C = 1); • nombre de paliers successifs sur lesquels la solution demeure inchangée: 5 (D = 1).

On suggère le centre du plan pour les paramètres A, B, C, et D, car se placer au niveau haut permet d'assurer un bon taux de convergence avec toutefois un nombre élevé d'évaluations de la fonction objectif (paliers de température longs, nombre élevé de paliers successifs pour l'arrêt de la procédure), alors que le niveau bas peut conduire à des taux de convergence assez faibles. Le centre du plan réalise donc un compromis acceptable.

  ∈ [-512, 512] f * -837. 966 (solution théorique) x1 * = x2 * = 420.97 (solutions théoriques) où les composantes x1 et x2 sont définies à 10 -3 près (Fig. 1).

Fig. 1 .- 3

 13 Fig. 1. Fonction BOITE.

Fig. 2 .

 2 Fig. 2. Fonction GRIE.

Fig. 3 .

 3 Fig. 3. Fonction Jong.

Fig. 4 .

 4 Fig. 4. Fonction Micha.

Fig. 5 .

 5 Fig. 5. Fonction Pib.

Fig. 6 .

 6 Fig. 6. Fonction Ros.

- 1 .

 1 25 -1.25 1.25 -3.75 1.25 -1.25 -1.25 -1.25 -1.25 -3.75 -3.75 1.25 3.75 3.75 1.25 Micha 1.25 -1.25 1.25 1.25 3.75 1.25 1.25 -6.25 8.75 -3.75 3.75 -1.25 1.25 3.75 -1.25 Pib -8.75 1.25 6.25 -3.75 1.25 6.25 -3.75 1.25 -3.75 6.25 -3.75 1.25 1.25 6.25 6.25 Ros -3.75 8.75 1.25 -3.75 -1.25 1.25 -3.75 -1.25 -1.25 -3.75 -1.25 -1.25 6.25 3.75 3.75

  10c) où n représente le nombre d'essais pour une combinaison donnée de paramètres (ici n = 5) et nev i est le nombre d'évaluations de la fonction objectif pour l'essai i. L'écart-type des estimations se déduit donc de celui des mesures (σ(sm2)) en le divisant par quatre. De fac ¸on classique, l'étude des variances et des intervalles de confiance est restreinte au centre du plan, qui correspond à A = 1, B = 1, C = 0.85 et D = 5.

  ). L'effet moyen le plus important constaté est dû au paramètre B, et le paramètre A a l'effet moyen le plus faible. En effet, le paramètre A représente la fac ¸on d'initialiser la température, et dans la mesure où le test d'arrêt est effectué sur la non évolution de la fonction objectif sur D paliers consécutifs et non sur une valeur minimale de la température,

	Tableau 11
	Valeurs des seconds members
	Boîte		Jong			Grie		Micha		Pib		Ros		Essai
	sm1	sm2	sm1	sm2		sm1	sm2	sm1		sm2	sm1	sm2	sm1	sm2
	0		100	3		80	6	20		0.4	100	7	0			1
	0		100	3		80	5	20		0.25	60	4	0			2
	40	5	100	36		80	54	20		0.3	100	78	20		14	3
	20	4.7	100	36		100	70	20		0.4	100	72	40		1.5	4
	0		100	14		100	31	60		2.8	8 0	3 0	2 0		7	5
	0		100	13		100	27	40		1	100	36	20		6	6
	60	26	100	198		100	358	0			100	390	40		7	7
	0		100	184		80	290	20		0.6	80	318	20		5.8	8
	20	15	100	6		100	1052	20		0.6	100	11	20		2	9
	0		6 0	4		8 0	9	2 0		0 .6	6 0	7	0			1 0
	20	6.6	100	45		80	89	40		4	80	77	40		1.7	1 1
	20	4.7	100	43		80	80	40		4	40	38	0			12
	40	3.8	100	29		80	126	20		0.2	100	63	0			13
	0		100	29		80	45	20		2.7	80	544	0			14
	40	3	100	226		80	350	20		14	100	460	20		85	15
	40	2.7	100	210		80	340	40		12	100	485	20		8	16
	Tableau 12
	Effets et interactions pour le critère sm1
	Fonction	A	B	C	D	AB	AC	AD	BC	BD	CD	ABC	ABD	ACD	BCD	ABCD
	Boîte	-8.75	11.25 3.75	3.75 -1.25 -3.75	1.25	1.25 -3.75	3.75 -1.25	8.75 1.25	1.25	3.75
	Jong	-2.5	2 .5	2.5	-2.5	2 .5	2 .5	-2.5	-2.5	2 .5	2 .5	-2.5	2 .5	2.5	-2.5	-2.5
	Grie

  des estimations, le nombre de degrés de liberté, la valeur de la variable de Student t, la valeur absolue h de l'intervalle de confiance pour sm2, et enfin la valeur absolue β de la borne de l'intervalle de confiance pour le rapport sm1/sm2. Les intervalles de confiance sont calculés au niveau de probabilité de 95%.Rappelons qu'un effet ou une interaction est significatif si son intervalle ne contient pas la valeur zéro; dans le cas contraire, cet élément est considéré comme non significatif. A partir des valeurs reportées dans le Tableau 13, sont mentionnés dans le Tableau 18 uniquement les éléments significatifs; les cellules vierges correspondent à des termes non significatifs.A titre d'exemple, l'intervalle de confiance pour l'effet de A pour la fonction Boîte qui est [-0.95 -1.17; -0.95 + 1.17] = [-2.12; + 0.22], permet de conclure que l'effet de A est non significatif.De même, l'intervalle de confiance pour l'effet de B pour la fonction Boîte, égal à [2.38 -1.17; 2.38 + 1.17] = [1.21; 3.55], indique qu'il est significatif.Cette étude statistique confirme et affine les résultats fournis par l'analyse qualitative, à savoir que l'effet de A est négligeable, celui de B est prédominant, C et D jouent des rôles comparables, avec un léger avantage pour C. Les interactions d'ordre deux faisant intervenir A ne sont pas significatives, ce qui est moins net pour les interactions d'ordre trois, et l'interaction d'ordre quatre est négligeable dans la majorité des cas.

	Tableau 15
	Valeurs de sm1 et sm2
	Fonction	sm1	sm2
	Boîte	40	3.182
	Jong	100	7.724
	Grie	80	13.578
	Micha	40	1.041
	Pib	90	42.973
	Ros	40	10.780
	des fonctions sont reportés pour chacun des 10 essais dans le
	Tableau 16, où les cellules vierges correspondent à des échecs
	(pas de convergence). Pour chaque fonction, le Tableau 17
	donne la moyenne sm2, l'écart-type σ(sm2), l'écart-type σ


Gwladys Ravon
email: gwladys.ravon@inria.fr

Florian Bernard

Angelo Iollo

Caroline Caradu

Automatic Rigid Registration of Aortic Aneurysm Arterial System

An abdominal aortic aneurysm is defined as a local and abnormal dilation of the aortic wall that can lead to rupture and death without treatment. A useful tool for the patient's postoperative follow-up is segmentation registration to see the evolution of the aneurysm between two examinations. Here, we propose a method to automatically register the entire arterial system: the segmentation is divided in three parts (suprarenal, infrarenal zone, and iliac arteries) and each part is registered separately. We chose a rigid point set to point set registration through the iterative closest point algorithm. We also compute the displacement fields and derive a criterion to accept or reject the registration of the infrarenal zone and iliac arteries. Registration is successful in 96% of cases for the infrarenal zone, in 94% for the suprarenal zone and in 65% for the iliac arteries.

Introduction

An abdominal aortic aneurysm is defined as a local and abnormal dilation of the aortic wall. Without treatment the aneurysm may grow until rupture, which ultimately leads to death in a majority of cases [START_REF] Chaikof | The Society for Vascular Surgery practice guidelines on the care of patients with an abdominal aortic aneurysm[END_REF][START_REF] Hoornweg | Meta analysis on mortality of ruptured abdominal aortic aneurysms[END_REF]. The preferred treatment modality now consists in the implantation of an endograft inside the aorta inserted through a sheath through the femoral arteries (endovascular aortic aneurysm repair (EVAR) [START_REF] Chaikof | Reporting standards for endovascular aortic aneurysm repair[END_REF]). The endograft isolates the aneurysm wall from the blood flow. Although it is a less invasive and morbid procedure procedure, the need for post-operative surveillance on a 1 long time scale is important [START_REF] Stather | Systematic review and meta-analysis of the early and late outcomes of open and endovascular repair of abdominal aortic aneurysm[END_REF]. The registration of information (images, centerlines or segmentations) obtained at different time points after the procedure is an important tool to follow the evolution of the aneurysm. It allows the clinician to see whether the aneurysm has grown or shrunk and where. Different approaches have been developed in the past years and we can cite the studies from Maiora et al. [START_REF] Maiora | Thrombus volume change visualization after endovascular abdominal aortic aneurysm repair[END_REF], Lopez-Linares et al. [START_REF] López-Linares | Image-based 3d characterization of abdominal aortic aneurysm deformation after endovascular aneurysm repair[END_REF], Braet et al. [START_REF] Braet | Vascular deformation mapping of abdominal aortic aneurysm[END_REF] or Demirci et al. [START_REF] Demirci | Quantification of abdominal aortic deformation after evar[END_REF]. These works have in common to combine several steps: they start with a rigid registration and the final step consists of a non-rigid deformation to align both objects and deduce the displacement field of the aneurysm between the two examinations. Another common point among these approaches is that at least one step of the process relies on the intensity of the images, thus the methods are not suitable for cases with no contrast agent. On the other hand, the basis for the rigid registration differ from one study to another: bones (Lopez-Linares et al.), lumen (Maiora et al.) or images (Braet et al. and Demirci et al.). The basis for the non rigid registration is also different: thrombus (Lopez-Linares et al.), lumen (Maiora et al.), images (Braet et al.) or centerlines (Demirci et al.). Methods are visually validated by experts, but a metric can be used to follow the improvement of the result at each step of the process. These four studies focus on the aneurysm in the infrarenal zone and neglect the suprarenal zone and the iliac arteries. In the present work, we propose a new method to register the arterial system based on a rigid point set registration. We chose a rigid transformation because we wanted to keep the aneurysm unchanged. Because of the patient's movement and breathing, the relative position of the arterial system may be different for different examinations, especially the outer parts. For this reason, in the following, segmentations were automatically divided into three parts (suprarenal, infrarenal zone and iliac arteries) and each part was registered independently. Cutting the segmentations was also useful to deal with different topologies: presence or absence of the aortic root, length of the iliac arteries, etc. Assuming a displacement in the normal direction we computed displacement fields and derived a criterion to accept or reject a registration without a visual validation.

Methods

Registration

Registration was performed based on an initial random distribution of points on the aortic wall (the outer surface of the segmentation) using the Iterative Closest Point (ICP) algorithm first introduced by Besl and McKay [START_REF] Besl | Method for registration of 3-D shapes[END_REF]. The method finds the transformation between two points clouds by minimizing the lowest distances between points (root mean square error). In the following, the number of points on each surface was set to 200. To improve the possibility of finding a global minimum we started by matching the centroids of the segmentations. In synthesis, by combing ICP and centroid matching, we determine an optimal rigid roto-translation that keeps the shape of the aneurysm unaltered.

Our objective is to study segmentations that are the result of an automatic workflow. Hence, as function of the image and the conditions of the scan setting, the data can correspond to different portions of the vascular network containing the aneurysm. Therefore, we initially performed the registration of the full segmentation which could go, for instance, from the lowest extremities to the aortic arch. However, the results were not satisfying since the ICP algorithm could not converge to a meaningful minimum if the branches of the two segmentations were too different. As a result, the registration of the full segmentation could not lead to an overall acceptable result.

The key point is to divide the segmentation into comparable regions. In this case, three parts: the suprarenal and infrarenal zones and the iliac arteries. Then, they were registered separately. As we wanted the workflow to be automatic, we developed our own cutting method (Riffaud et al, submitted). This method identifies the different branches of the arterial system from the abdominal aortic aneurysm segmentation, through the centerline and its bifurcation points. The application of this method to our segmentation identified the bifurcations of the renal arteries and the aortic bifurcation. Once the vascular branches are correctly identified, the segmentation is cut above the highest renal artery and 30 mm below the aortic bifurcation. Retention of the renal and incipient iliac arteries was suitable to ensure the finding of the correct rotation.

Additional preprocessing was required for the three regions. First, we needed to remove outliers from the renal arteries for the infrarenal zone. Regarding the suprarenal zone, both the celiac artery and the superior mesenteric arteries led the ICP algorithm to a local minimum rather than the global minimum. It was even less effective in presence of complex curved vessels. We had to remove them to ensure a good result. To this end, we applied an erode filter to both segmentations. To avoid problems caused by the presence of the aortic arch in only one segmentation, we kept the same height of the aorta for the target and the source. It was determined as the smallest height between the 2 suprarenal zones. We performed the registration of the two corresponding surfaces and applied the transformed matrix to the original source surface.

Regarding the iliac arteries, we registered separately the right and left common iliac arteries to minimize the impact of patient's motion. To ensure convergence to a global minimum we kept the same length on both segmen-

Displacement fields

We computed the displacement fields in the normal direction of the aneurysm between two scans. For each point on the target surface, we determined if it was inside (respectively outside) the registered source surface. We then found the closest intersection point with the registered source in the outward (resp. inward) normal direction, at most at 20 millimeters. The positive (resp. negative) euclidean distance between the two points determined the value of the displacement. An example can be seen in Figure 2: we represented displacement fields from dark blue if the aneurysm had significantly

Accepted and rejected registrations

All the registrations were visually evaluated by experts and registration was rejected when one could notice a left-right or top-bottom shift, or when it was visible that the rotation was wrong. In the figures white surfaces are the target surfaces, obtained from the first post-EVAR examination. Red surfaces are the registered source surfaces from the second examination. Figure 3 shows examples of rejected registrations of the full geometry: for Patients A and C there were lateral shifts, whereas for Patient B it was a vertical shift. We show examples of accepted registrations of the infrarenal zone for the three same patients in Figure 4. Figure 5 exhibits rejected registrations of suprarenal zones for the same patients: mainly the rotation was not good. Figure 6 shows accepted registration of the same patients with our proposed method. Concerning the iliac arteries Figure 7 shows accepted registrations for Patient A and B even though the end of the iliac arteries did not match well. This was expected due to the changes in the patient's position. For Patient C, the registration was rejected because of a slight sub-shift.

Criterion

With the goal of automatically distinguishing an accepted registration from a rejected one, we searched for a criterion based on the statistics of the values of the displacement field. To evaluate the registration of the full geometry and the infrarenal zone we automatically select a centered 10-millimeter-thick slab at the point with the largest axial diameter (see Figure 8). For the suprarenal zone, we considered the lowest 30 millimeters of the segmentation and the 30 highest for the iliac arteries. We then computed the median (M ) and the first (Q 1 ) and third (Q 3 ) quartiles of the displacement field values on these slabs. We introduce the product

P = Q 1 * Q 3 . (1)
This product P indicates if there is a change in the sign of the displacement field, which could correspond to a shift in the registration. On the other hand, if the product is positive it means that the aneurysm has evolved in the same way on the full slab: either shrunk or grown. We also define the dispersion D as

D = Q 3 -M. (2)
This quantity gives an indication of the dispersion of the displacement field values. In the next section, the objective is to find a threshold value for D to accept or reject a registration, depending on the sign of P.

Dataset

The dataset was composed of segmentations from 42 different patients for a total of 81 registrations with always the segmentation obtained right after the procedure as the target segmentation. All the segmentations were provided by Nurea (https://www.nurea-soft.com). The methodology was validated by clinicians and was the object of the study published by Caradu et al. [START_REF] Caradu | Fully automatic volume segmentation of infrarenal abdominal aortic aneurysm computed tomography images with deep learning approaches versus physician controlled manual segmentation[END_REF].

Results

All the registrations were visually labeled as accepted or rejected by experts.

Full geometry and infrarenal zone

Among the 81 registrations of the full geometries, only 11 were accepted. Among the 70 rejected 36 were of really poor quality with no superposition of the target and the registered source (see Figure 9 for an instance). Most of the time this was due to an important difference in the topology like the presence or absence of the aortic arch. When we computed the displacement fields, it was initialized to zero, so if the registered source was too far from Regarding the registrations of the infrarenal zones, 78 were accepted (96%) and none were of really poor quality. There was a problem with the segmentation in the 3 rejected cases: artifacts or superior mesenteric artery glued to the sac. Figure 10 shows the distribution of D depending on the sign of P for rejected and accepted registrations. One can see that a threshold can be found to distinguish rejected registrations from accepted ones, especially when the product P is negative. To find those thresholds we used a Support Vector Machine (SVM) algorithm with a linear kernel on two sets: positive product on one hand, negative product on the other hand. We then predicted the output (accepted or rejected) for a range of dispersion values between 0.1 and 10 with a step of 0.1. The results can be seen in Figure 11: the threshold for a negative product was much smaller than for a positive product, respectively 1.6 and 3.54. If we applied these thresholds on our registrations, 11 would have been wrongly accepted (10 full geometries and 1 infrarenal zone) and 2 would have been wrongly rejected (2 infrarenal zone).

In short, a registration should be rejected in the following cases:

• P = 0,

• P < 0 and D > 1.6,

Suprarenal zone

For the suprarenal zone, only 5 registrations were rejected including 3 with artifacts and 2 with a thin aorta. Only 3 cases had a positive product P so we did not make any distinction according to its sign. In Figure 12 we can see that no distinction can be made between values for accepted and rejected registrations. When we applied an SVM algorithm, the threshold between accepted and rejected was 86 in the linear case and 135 in the nonlinear case which was not relevant. Moreover, when we removed the isolated accepted case, the predictions always resulted in "accepted" in both cases, which may be due to the small number of rejected cases. We can assume that the registration is accepted if the segmentations are correct.

Iliac arteries

Registration of this part was the most difficult because of the shape of the iliac arteries: 53 registrations were accepted (65%). 20 were rejected with only one common artery correctly registered and 8 were rejected with both arteries badly registered. Among those 28 cases, 7 presented artifacts in the segmentation or had at least one very short common iliac artery. The success rate rises to 72% when the segmentations are clean. Only 6 cases had a positive product so it was meaningless to make a difference. The SVM algorithm on the distribution of the dispersion D visible in Figure 13 Figure 12: Dispersion D for registrations of the suprarenal zones won't be relevant. However, a dispersion D greater than 7 mm indicates a non-satisfying registration. 9 cases had a null product P, including 4 with artifacts or short iliac arteries. In short, a registration of the iliac arteries should be rejected in the following cases:

• P = 0,

• P = 0 and D > 7.

13 registrations would have been rightly rejected while 4 would have been wrongly rejected. In 2 cases the source was cut lower than the target so that only a small portion of the registered source was within the 30 mm considered for the computation of the dispersion D.

Figure 14 presents the results with the three different parts on the same image for Patients A, B and C. Small discontinuities can be seen between parts on Patients B and C.

Evolution in time: example of patients' follow-up

Follow-up of the surgery includes regular examinations to control the evolution of the aneurysm. That is why we had multiple segmentations for a given patient and could perform registration by parts of these segmentations. The target surface was always the surface from the scan performed immediately after the procedure. For Patient A we had two segmentations obtained after 3 months and 1 year. In Figure 15 one can see that the aneurysm started to shrink after 3 months and the reduction was more important after one year. The aneurysm was also partially located in the right common iliac artery, For Patient B we had 4 scans: 3 and 6 months and 1 and 2 years after the procedure. Figure 16 shows that the aneurysm remained stable after 3 months. It then started to grow at 6 months, remained stable between 6 month and 1 year but then had a massive global growth after 2 years. With the displacement field (Figure17) we retrieved the stable-growth-stable-massive growth scheme.

For Patient C we had only one scan 4 months post-EVAR and the result is visible in Figure 14. The aneurysm seemed to be stable. However with the visualization of the displacement field (Figure 18) we can see that the aneurysm had grown. It shows that the aneurysm had grown locally on the back side of the sac while the front side remained stable. In this particular case, the maximum diameter had not changed despite an increasing volume. We have developed an automatic pipeline that registers the arterial system and allows visualization of an aortic aneurysm evolution for a patient between two examinations. Unlike other studies, we used a rigid registration and assumed that the aneurysm evolves in the normal direction to the wall of the aneurysm.

Registering the entire arterial system up to the iliac arteries, rather than just the aneurysm region, makes it possible to follow any aneurysm regardless of its location. We proposed a method based on the iterative closest point algorithm to register segmentations of abdominal aortic aneurysms. It registers not only the aneurysm region but the entire arterial system including the aorta and iliac arteries. Segmentations were divided into 3 parts: suprarenal, infrarenal zones and the iliac arteries. Each region was registered separately, after a preprocessing depending on the region. We have shown that registration of the infrarenal zone instead of the whole segmentation leads to better results in terms of expert assessment. Our method leads to accepted registrations in 96%, 94% and 65% of cases for the infrarenal and suprarenal zones and the iliac arteries, respectively. 78% of individual iliac arteries are registered in an acceptable manner.

Rigid registration allows an objective, model-free, comparison to evaluate the evolution of aneurysm surfaces. Clearly, the limitations of our approach are intrinsic to the minimization algorithm. In extreme situations, for example in the case of large variations in tortuosity, minimization will somehow average the displacement of the centerline and that of the surfaces. We have also introduced the computation of displacement fields in the normal direction. The displacement field represents a quantitative indicator for the evaluation of the aneurysm's evolution. Its visualization is quite simple with the color scale. Based on its value, we have also developed a criterion to accept or reject the registration of the infrarenal part and the iliac arteries. None was required for the supra-renal zone. Our method, unlike other existing methods, is independent of the images or anatomic location of the aneurysm. It can be performed with or without contrast product and is not disturbed by the presence of artifacts in the image. It can also be interpreted as a preliminary step that can be used as a first step for non-rigid registration.

The choice of a point set registration algorithm was guided by the imperative of user-free usage. For example, a method that relies on landmark registration usually requires user intervention to select or validate landmarks which also involve uncertainty and inter-user variability. Moreover, the fewer landmark points considered, the more important the accuracy of the position of each point. In contrast, this problem is avoided with the aortic wall surfaces due to the systematic random sampling of the surfaces.

A known bias in aneurysm volume measurement is the cardiac cycle [START_REF] Truijers | In-vivo imaging of changes in abdominal aortic aneurysm thrombus volume during the cardiac cycle[END_REF][START_REF] Grøndal | The cardiac cycle is a major contributor to variability in size measurements of abdominal aortic aneurysms by ultrasound[END_REF]: it is not constant during the whole cycle. Part of the changes in volume or shape we observed could be attributable to a different phase of the cardiac cycle. However, this is impossible to quantify without having full control on the image acquisition. To avoid this bias, it would be necessary to register segmentations taken during the same phase. However, such information is not available from CT-scan images.

Registration is an important tool to follow the evolution in time of an aneurysm. As an example of application, we showed the results for 3 patients: one with a shrinking aneurysm (Patient A), one with a globally growing aneurysm (Patient B) and one with a locally growing aneurysm (Patient C). For Patient A, the aneurysm was also located in one common iliac artery. Thanks to our method we could follow its evolution and see that it had shrunk in the iliac artery too.

Figure 1 :

 1 Figure 1: Schematic representation of the automatic registration pipeline

Figure 2 :

 2 Figure 2: Displacement fields for Patient A 3 months and 1 year after the procedure. The target surface was the one from the post-operative scan

Figure 3 :

 3 Figure 3: Example of registration of the full surface for 3 patients. White surfaces are the target surfaces, obtained from the first post-EVAR examination. Red surfaces are the registered source surfaces from the second examination

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: Registration of the infrarenal zone for Patients A, B and C. Registered infrarenal zone is superimposed on the full geometry

Figure 8 :

 8 Figure 8: Example of a slab around the point with the highest axial diameter. Segmentation (blue), centerline colored with the diameter values and 10 millimeters thick slab

Figure 9 :

 9 Figure 9: Example of rejected and very poor registration of full segmentations

Figure 10 :Figure 11 :

 1011 Figure 10: Dispersion D for registrations of the full geometries and the infrarenal zones. Top: P > 0, bottom: P < 0. Red: rejected registrations; blue: accepted ones

Figure 13 :

 13 Figure 13: Dispersion D for registrations of the iliac arteries

Figure 14 :

 14 Figure 14: Registration of the three different parts for Patients A, B and C

Figure 16 :

 16 Figure 16: Registrations for Patient B for 4 scans after the procedure. The target surface was the one from the post-operative scan

Figure 17 :

 17 Figure 17: Displacement fields for Patient B for 4 scans after the procedure. The target surface was the one from the post-operative scan

Author contribution

GR designed the study, implemented the algorithms, analyzed and interpreted the results and drafted the manuscript. FB and AI helped conceive the study, provided feedback about the results and revised the manuscript. CC provided expert analysis on the results and revised the manuscript.

Conflict of interest

The authors declare that they have no conflict of interest.

G. Ravon completed a PhD in applied mathematics at the University of Bordeaux. The subject mixed inverse problems and cardiac electrophysiology. She has currently a post-doctoral position at Inria Bordeaux Sud-Ouest in Memphis Team. Her main research topic is applied mathematics to medicine and health care.

F. Bernard completed a PhD in applied mathematics at University of Bordeaux and a PhD in fluid dynamics at Politecnico di Torino in 2015. He's now CEO of Nurea, a company dedicated to the development of decision making support software for physicians to improve diagnostic, follow up and pronostic of cardiovascular diseases.

A. Iollo is professor of Applied Mathematics at the University of Bordeaux and head of team Memphis at Inria, the national institute for applied mathematics and computer science.

Dr C. Caradu is a M.D. Ph.D. working as a vascular surgeon at the University Hospital of Bordeaux. Her research interests include the endovascular treatment of complex aortic aneurysms, the use of artificial intelligence in the treatment of vascular pathologies, aortic infections, and the role of endothelial cell dysfunction in the pathophysiology of critical limb ischemia.