--- title: VISI601_CMI TP3 Calcul de distances et de plus courts chemins via des équations différentielles --- # TP3 - Calcul de distances et de plus courts chemins via des équations différentielles > [name=Jacques-Olivier Lachaud][time=February 2022][color=#907bf7][time=March 2026][color=#907bf7] > (Les images peuvent être soumises à des droits d'auteur. Elles sont utilisées ici exclusivement dans un but pédagogique) ###### tags: `visi601` `tp` Retour à [VISI601_CMI (Main) Algorithmique numérique](https://codimd.math.cnrs.fr/s/IWTaBkA9m) --- [TOC] --- :::warning Archive d'images à récupérer : [TP3.zip](https://jacquesolivierlachaud.github.io/teaching/algo-num/Tests/TP3.zip) [TP3.tgz](https://jacquesolivierlachaud.github.io/teaching/algo-num/Tests/TP3.tgz) ::: --- ## Introduction Ce TP va vous montrer comment calculer les distances à des sources placées dans une grille avec des obstacles, ainsi que des plus courts chemins. Notre approche va utiliser des équations différentielles simples, basées sur l'opérateur Laplacien. Elle est une adaptation de la méthode *[Geodesics in heat](https://jacquesolivierlachaud.github.io/lectures/info601cmi/Cours/crane-2013.pdf)* de Crane, Weischedel et Wardetsky. Comme nous allons nous placer dans le plan, l'opérateur Laplacien sera plus simple à écrire que dans l'article. Sur le fond, nous allons faire la même chose. Dans un domaine $\Omega$, on cherche à calculer la distance (et aussi la direction du plus court chemin) à un ensemble $X$ donné. L'idée est de diffuser de la chaleur à partir de $X$, puis de renormaliser les lignes équidistantes de la chaleur. 1. on définit $X \subset \Omega$ comme des sources de chaleurs, donc $u_0(x,y)=1$ si $(x,y) \in X$ et $u_0(x,y)=0$ sinon. 2. on diffuse la chaleur pour un temps $T$ avec l'équation de la chaleur $\partial_t u = \Delta u$, en partant de $u(0)=u_0$ 3. on regarde les directions de variation de la chaleur $u$ au temps $T$. On calcule donc le champ de vecteurs $V = − \frac{\nabla u}{|\nabla u|}$. Chacun de ses vecteurs (tous de norme 1) pointe dans la direction où on s'éloigne le plus vite de la chaleur (donc de la source de chaleur la plus proche). 4. enfin on résoud l'équation de Poisson $\Delta d = \mathrm{div} V$, où $\mathrm{div} V$ est la divergence de $V$, égale à $\frac{\partial V_x}{\partial x}+\frac{\partial V_y}{\partial y}$. |  | | ----------------------------------------------------------------------------------------------------- | | Figure 12 de l'article: en haut on voit l'influence du choix des conditions aux bords pour le Laplacien dans les plus courts chemins (gauche : Neumann, droite : Dirichlet), en bas, on voit l'influence du temps de diffusion (petit $T$ à gauche, grand $T$ à droite) | Nous allons mettre en oeuvre cette méthode dans un domaine du plan défini par une image. On pourra utiliser les images suivantes, où le domaine $\Omega$ correspond aux pixels de l'image qui ne sont pas noirs. |  |  |  | | --------------------------------------------------------------------------------------------------- | --------------------------------------------------------------------------------------------------- | --------------------------------------------------------------------------------------------------- | > Un informaticien pourrait se demander pourquoi on n'utilise pas un algorithme de calcul de plus courts chemins dans un graphe, comme l'algorithme de Dijkstra. En fait, si on définit le graphe avec les 4 voisins de chaque pixels, ces algorithmes retournent le plus court chemin selon la distance de Manhattan (comme si on se déplace sur un échiquier). Le résultat ne serait donc pas du tout précis ou isotrope. [color=#907bf7] > Dans tout le TP, on caractérise un pixel par ses coordonnées $(i,j)$, où $i$ est sa ligne et $j$ sa colonne. C'est donc la convention usuelle des matrices, mais pas celle des coordonnées dans le plan, qui est plutôt $(x,y)$ où $x$ correspond à l'abscisse (donc une colonne) et $y$ à l'ordonnée (donc une ligne). [color=#cc0000] **A rendre** : Vous me ferez une archive contenant le ou les programmes python, ainsi qu'un compte-rendu donnant l'état d'avancement et les différentes choses demandées dans l'énoncé. Il sera tenu compte dans la notation du respect des consignes, de la qualité de code, mais aussi de la façon dont chaque partie du travail demandé est testée (dans le compte-rendu ou dans le code lui-même). Vous me remettrez l'archive via <a href="https://tplab.apps.math.cnrs.fr">TPLab</a>, en fin de séance, puis une version finale avant le **dimanche 26 avril, minuit**. ## I Définition du Laplacien pour un domaine quelconque Vous allez reprendre votre classe `Grid` définie dans le TP2 précedent. On va faire un nouveau constructeur de cette classe, qui ne va définir des indices que sur les pixels qui ne sont pas noirs (i.e. non nuls). On rappelle que pour charger une des images, on peut écrire: ```python= # load color image img=mpimg.imread('exemple-1.png') red = img[:, :, 0] # image N&B entre 0 and 1, qui correspond au canal rouge imgplot = plt.imshow(red,cmap='gray',vmin=0.,vmax=1.) plt.show() # éventuellement ``` La variable `red` est un `numpy.array`, qui correspond à une matrice de la hauteur de l'image et de la largeur de l'image. Le nouveau constructeur va donc parcourir l'image pour trouver les pixels "utiles" (i.e. ceux qui appartiennent au domaine $\Omega$), les numéroter, et mettre leurs coordonnées (i,j) dans les données membres `self.I` et `self.J`. Il faut aussi remplir le dictionnaire `index` de manière à relier des coordonnées `(i,j)` à un numéro de pixel `idx`. :::info **Complétez** donc le constructeur ci-dessous. ```python= import numpy as np from scipy.sparse import * import matplotlib.pyplot as plt import matplotlib.image as mpimg class Grid: ... def __init__( self, image ): self.nbrow = image.shape[0] self.nbcol = image.shape[1] self.I = [] self.J = [] self.index={} # à compléter ``` ::: Maintenant, on crée une grille spécifique à l'image chargée ainsi: ```python= G=Grid(red) ``` Vérifiez maintenant que vous avez bien défini votre Laplacien, en regardant pour chacun des 4-voisins si il existait via `Grid.getIndex(i+1,j)`, etc. Normalement vous n'avez aucun code à modifier !! Si vous testez ainsi ```python= img=mpimg.imread('exemple-1.png') red = img[:, :, 0] # image N&B entre 0 and 1, qui correspond au canal rouge G=Grid(red) L=G.Laplacian() # Laplacien associé à Neumann L ``` doit afficher la matrice creuse ainsi: ```shell <114842x114842 sparse matrix of type '<class 'numpy.float64'>' with 568438 stored elements in Compressed Sparse Column format> ``` Vous pouvez vérifier votre code sur le petit exemple ci-dessous: ```python= A=np.zeros( (5,5)) A[:,1]=1. A[2,:]=1. A[3,2]=1. print(A) G=Grid(A) L=G.Laplacian() # Laplacien pour les conditions de Neumann print(L) ``` qui donne ```shell [[0. 1. 0. 0. 0.] [0. 1. 0. 0. 0.] [1. 1. 1. 1. 1.] [0. 1. 1. 0. 0.] [0. 1. 0. 0. 0.]] (0, 0) -1.0 (1, 0) 1.0 (0, 1) 1.0 (1, 1) -2.0 (3, 1) 1.0 (2, 2) -1.0 (3, 2) 1.0 (1, 3) 1.0 (2, 3) 1.0 (3, 3) -4.0 (4, 3) 1.0 (7, 3) 1.0 (3, 4) 1.0 (4, 4) -3.0 (5, 4) 1.0 (8, 4) 1.0 (4, 5) 1.0 (5, 5) -2.0 (6, 5) 1.0 (5, 6) 1.0 (6, 6) -1.0 (3, 7) 1.0 (7, 7) -3.0 (8, 7) 1.0 (9, 7) 1.0 (4, 8) 1.0 (7, 8) 1.0 (8, 8) -2.0 (7, 9) 1.0 (9, 9) -1.0 ``` ## II Initialisation des sources On va créer deux méthodes pour initialiser le vecteur $\mathbf{u}_0$, qui doit valoir 1 sur les pixels sources du domaine (i.e. les points de $X \subset \Omega$). 1. Une méthode `Grid.sources` qui fabrique le vecteur `U0` à partir de la liste des coordonnées des pixels sources 2. Une méthode `Grid.boundary` qui fabrique le vecteur `U0` en considérant que les sources sont les pixels du bord du domaine (i.e. ceux qui n'ont pas 4 voisins dans le domaine). :::info Ecrivez ces 2 méthodes: ```python= class Grid: ... # given a list of pixels [(i,j),...], return the vector U0 with a value 1 on these sources def sources( self, H ): ... # return U0 where its value is 1 on its boundary pixels def boundary( self ): ... ``` ::: On peut maintenant définir l'initialisation pour une source en (1,1): ```python= U0=G.sources([(1,1)]) ``` ou comme des sources sur tous les bords du domaine: ```python= U0 = G.boundary() Visu_U0 = G.vectorToImage( U0 ) imgplot = plt.imshow(Visu_U0,cmap='gray',vmin=0.,vmax=1.0) ```  ## III Diffusion de la chaleur à partir des sources On va maintenant diffuser le vecteur $\mathbf{u}_0$ comme si c'était des sources de chaleur qui se propagent. Il suffit donc d'utiliser une méthode de diffusion du tp précédent pour calculer cette diffusion. La méthode d'Euler implicite permet de faire des pas beaucoup plus grands. :::info Vérifiez que votre diffusion implicite fonctionne bien avec: ```python= # diffuse jusqu'à T=1000, par pas dt=100 U=G.implicitEuler( U0, 1000, 100 ) # T=1000 visu_U = G.vectorToImage( W ) # notez le rescale à 1e-6 pour bien voir la "chaleur" plt.imshow(visu_U,cmap='hot',vmin=0,vmax=1e-6) ``` ::: On voit que le temps de diffusion influe sur la propagation de la chaleur à partir de la source donnée. |  |  |  | | -------- | -------- | -------- | | $\mathbf{u}_T$, $T=100$ | $\mathbf{u}_T$, $T=1000$ | $\mathbf{u}_T$, $T=10000$ | ## IV Calcul du gradient normalisé De la même façon qu'on a définit un opérateur Laplacien, nous allons définir deux opérateurs pour le gradient $\nabla u = \left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}\right)$, l'un pour calculer la dérivée selon $x$, l'autre selon $y$. On va donc fabriquer deux matrices creuses $G_x$ et $G_y$ qui approchent les deux composantes du gradient. Ainsi $G_x \mathbf{u}$ sera le vecteur stockant la composante $x$ du gradient pour tous les pixels du domaine, et $G_y \mathbf{u}$ la composante $y$. Pour un pixel $(i,j)$, il peut y avoir 4 cas pour approcher le gradient selon $x$ en ce point: 1. si $(i,j-1) \in \Omega$ et $(i,j+1) \in \Omega$ alors $\frac{\partial u}{\partial x}(i,j) \approx (\mathbf{u}(i,j+1) - \mathbf{u}(i,j-1))/2$ 2. si $(i,j-1) \not\in \Omega$ et $(i,j+1) \in \Omega$ alors $\frac{\partial u}{\partial x}(i,j) \approx \mathbf{u}(i,j+1) - \mathbf{u}(i,j)$ 3. si $(i,j-1) \in \Omega$ et $(i,j+1) \not\in \Omega$ alors $\frac{\partial u}{\partial x}(i,j) \approx \mathbf{u}(i,j) - \mathbf{u}(i,j-1)$ 4. et sinon $\frac{\partial u}{\partial x}(i,j) \approx 0$. Pour approcher le gradient selon $y$, le principe est similaire, mais selon les lignes. :::info Ecrivez donc les méthodes de `Grid` qui calculent les deux opérateurs gradient $G_x$ et $G_y$. ```python= class Grid: ... def gradX(self): ... def gradY(self): ... ``` ::: :::info Ensuite utiliser ces 2 opérateurs pour calculer l'opposé du gradient normalisé, i.e. $− \frac{\nabla u}{|\nabla u|}$. La méthode `direction` doit donc retourner une paire de vecteurs colonne `Vx,Vy` à partir du vecteur colonne `U`. Il faudra donc diviser les vecteurs par leur norme pour obtenir des vecteurs unitaires. ```python= class Grid: ... def direction(self,U): ... ``` ::: > On rappelle que `numpy` permet de faire des opérations sur tous les éléments d'un tableau `X` comme `np.pow(X,2.0)`... Vous pouvez vérifier votre méthode sur l'image `exemple-1.png`, $T=100$. ```python= VX,VY = G.direction( U ) visu_VX = G.vectorToImage( VX ) plt.imshow(visu_VX,cmap='viridis',vmin=-1.,vmax=1.) visu_VY = G.vectorToImage( VY ) plt.imshow(visu_VY,cmap='viridis',vmin=-1.,vmax=1.) ``` |  |  | | -------- | -------- | | Composante $x$ de $(V_x,V_y)$ | Composante $y$ de $(V_x,V_y)$ | ## IV Résolution de l'équation de Poisson Nous avons donc déterminé un **champ de vecteurs** $V=(V_x,V_y)$ qui associe à tout pixel un vecteur unitaire donnant la direction où on s'éloigne le plus vite des sources. *Un déplacement de longueur 1 dans cette direction augmente donc la distance de 1 par rapport aux sources*. :::warning Idéalement, si $d$ est une fonction donnant la distance aux sources, on devrait avoir la relation: $\forall (x,y) \in \Omega, \nabla d (x,y) =V(x,y)$, c'est-à-dire que le gradient de la distance est exactement notre champ de vecteurs $V$. En pratique notre champ de vecteurs n'est pas exact, et n'est probablement pas le gradient d'une fonction. ::: En conséquence, plutôt que d'essayer résoudre $\nabla d = V$ (qui conduit à des systèmes sur-contraints sans solution exacte), on va résoudre un problème plus simple qui lui aura des solutions. On prend donc la divergence sur les 2 membres de l'équation précédente, et on aboutit à: $\forall (x,y) \in \Omega, \mathrm{div} \nabla d(x,y) = \mathrm{div} V(x,y) \quad \Longleftrightarrow \quad \forall (x,y) \in \Omega, \Delta d(x,y) = \mathrm{div} V(x,y)$. On rajoute la condition aux bords que $\forall (x,y) \in X, d(x,y)=0$. Si on discrétise l'équation précédente, on retrouve notre matrice Laplacien $L$ à gauche, et il nous reste juste à calculer le membre droit. :::info La **divergence** du champ $V$ s'écrit simplement car il s'agit de la somme du gradient selon $x$ de $V_x$ et du gradient selon $y$ de $V_y$. Comme vous avez déjà écrit les fonctions `Grid.gradX` et `Grid.gradY`, la méthode `Grid.divergence` s'écrit en 1 ligne ! ```python= class Grid: ... def divergence( self, Vx, Vy ): ... ``` ::: Soit $\mathbf{b} = \mathrm{div} V$, il ne vous reste plus qu'à résoudre l'**équation dite de Poisson** $L \mathbf{d} = \mathbf{b}$ où $L$ est votre opérateur matriciel approchant le Laplacien. :::info Ecrivez la méthode `Grid.poisson` qui résoud l'équation précédente. On peut utiliser la méthode PLU, par exemple via `scipy.linalg.splu`. ```python= class Grid: ... # Résoud Lx = B def poisson( self, B ): ... ``` ::: > Attention, le laplacien $L$, quoiqu'une matrice $n \times n$, n'est pas de rang maximum (car le vecteur constant $\mathbf{1}$ induit $L \mathbf{1} = \mathbf{0}$), et donc non inversible. Celà ne gêne pas PLU, mais la solution est correcte à une constante près quelconque. Il suffit donc de soustraire le minimum de la solution $\mathbf{d}'$ trouvée à chaque élement de $\mathbf{d}'$ pour trouver la solution $\mathbf{d}$, qui vaut bien 0 sur les sources. [color=#cc0000] On obtient les résultats suivants, pour la source en (1,1) et différents temps $T$. ```python= B = G.div( UX, UY ) D = G.poisson( B ) visu_D = G.vectorToImage( D ) plt.imshow(visu_D,cmap='hot') ``` |  |  | | | -------- | -------- | -------- | | $\mathbf{d}$, $T=100$ | $\mathbf{d}$, $T=1000$ | $\mathbf{d}$, $T=10000$ | :::info Calculez les distances sur d'autres exemples. Est-ce que vous constatez une influence du temps $T$ ? ::: ## V Influence de l'opérateur Laplacien On voit que l'opérateur Laplacien intervient deux fois dans la méthode 1. la première fois lors de la diffusion à partir des sources, 2. la deuxième fois lors de la reconstruction de la distance par la résolution de l'équation de Poisson. On a définit deux opérateurs Laplacien, $L_N$ imposant les conditions dites de Neumann (`Grid.Laplacian`) et $L_D$ imposant les conditions dites de Dirichlet (`Grid.LaplacianD`). Pour vérifier plus précisément l'influence de ces choix, je vous donne une fonction qui permet d'afficher la fonction distance de manière plus jolie. On voit ainsi mieux les différences entre calculs de distance selon $T$, ou selon l'opérateur Laplacien. ```python= # n est le nombre de lignes de niveaux visibles, width est leur largeur. def transformDistance( D, n = 50, width = 0.2 ): TD = D.copy() M = D.max() E = M / n for i in range( D.shape[0]): v = D[i] / E e = v-floor(v) if ( e < width ): TD[ i ] = 1.5*M return TD ``` :::info **Vérifiez** que ces opérateurs influencent grandement les résultats. Voilà ce que vous devez obtenir ($T=100$, $dt=100$, 15 isolignes de 0.2, image `exemple-4.png`, source en $(200,200)$): | $\mathbf{u}_T$, $T=1000$ | Diffusion par $L_N$ | Diffusion par $L_D$ | | -------- | -------- | -------- | | Poisson par $L_N$ |  | | | Poisson par $L_D$ |  |  | ::: Les approches où on ne prend pas le même Laplacien dans les deux opérations ne donnent pas forcément des résultats très consistants. :::warning Si on prend $L_N$ (Neumann) pour diffusion et reconstruction, alors les isolignes de distance ont tendance à être **perpendiculaires** aux bords. Si on prend $L_D$ (Dirichlet) pour diffusion et reconstruction, alors les isolignes de distance on tendance à être **parallèles** aux bords. ::: :::success Les auteurs proposent donc de prendre la **moyenne** des deux résultats pour **ignorer** les effets de bords sur la reconstruction de la distance. On va donc trouver $\mathbf{d}_N$ la distance en utilisant l'opérateur $L_N$ pour diffusion et Poisson, et on va trouver $\mathbf{d}_D$ la distance en utilisant l'opérateur $L_D$ pour diffusion et Poisson, puis calculer leur moyenne $\frac{\mathbf{d}_N + \mathbf{d}_D}{2}$. --- $T=100$, $dt=100$, 15 isolignes d'épaisseur 0.2 | $\mathbf{d}_N$ | $\mathbf{d}_D$ | $\frac{\mathbf{d}_N + \mathbf{d}_D}{2}$ | | ----------------------------------------------------------------------------------------------------------------- | ----------------------------------------------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------- | |  |  |  | |  |  |  | ::: ## VI Méthode "Heat" au complet On peut mettre maintenant toutes ces fonctions successives ensembles pour calculer en un appel la fonction distance à des sources. :::info Ecrivez donc la méthode `Grid.heatMethod` qui prend en entrée le vecteur `U0`, le laplacien choisi ('Neumann' ou 'Dirichlet'), le temps `T` et le pas `dt`. ```python= class Grid: ... def heatMethod( self, U0, Lap='Neumann', T=100, dt=10 ): if ( Lap == 'Neumann' ): self.L = self.Laplacian() else: self.L = self.LaplacianD() # à compléter ``` On note que ça vaut le coup de mémoriser une fois pour toutes les opérateurs Laplacien $L$ et gradients $G_x$ et $G_y$ car ils servent à plusieurs endroits. ::: On peut maintenant afficher la distance `D` calculée transformée pour voir les lignes de niveaux de cette fonction. ```python= TD = transformDistance(D,50,0.2) visu_TD = G.vectorToImage( TD ) plt.imshow(visu_TD,cmap='hot') ``` --- Solution pour le **Laplacien avec conditions de Neumann**. | |  |  | | -------- | -------- | -------- | | $\mathbf{d}_N$, $T=100$ | $\mathbf{d}_N$, $T=1000$ | $\mathbf{d}_N$, $T=10000$ | --- Si on suit les recommandations des auteurs, on peut faire la moyenne des solutions pour le **Laplacien avec conditions de Neumann** et **Laplacien avec conditions de Dirichlet**. On voit qu'un petit $T$ donne un résultat plus précis de la fonction distance, un plus grand $T$ fabrique des distances plus lisses. En revanche, si on descend trop $T$, on induit des divisions par zéro pour des diffusions trop courtes. |  |  |  | | -------- | -------- | -------- | | $\frac{1}{2}(\mathbf{d}_N+\mathbf{d}_D)$, $T=100$ | $\frac{1}{2}(\mathbf{d}_N+\mathbf{d}_D)$, $T=1000$ | $\frac{1}{2}(\mathbf{d}_N+\mathbf{d}_D)$, $T=10000$ | --- :::info Calculez les transformées en distance sur les autres images proposées. Essayez de mettre d'autres sources pour valider votre méthode. ::: ## VI Squelette ou axe médian d'une forme L'axe médian est le lieu des points d'un ensemble $\Omega$ tel que ce point soit à égale distance d'au moins 2 points du bords. Ca correspond aussi à des discontinuités de la fonction distance, ou des lignes de crête sur cette fonction distance (vue comme carte d'élévation). Ca ressemble aussi à un squelette de la forme, et cela le rend très utile en analyse de formes. Ce n'est pas si facile que ça de la calculer, mais on se propose ici de l'approcher en extrayant les lignes de crête de la fonction distance. Nous suivons l'approche de [Ridge Detection (Wikipedia)](https://en.wikipedia.org/wiki/Ridge_detection). :::info Vous définirez une méthode `Grid.ridges( self, D )` qui, étant donné le vecteur des distances $D$ aux bords du domaine, retourne un vecteur de valeur valant 1 sur les lignes de crête, et 0 sinon. On suivra les étapes suivantes 1. il faut lisser un peu la fonction distance $D$ (que vous avez calculée ci-dessus) pour que ses dérivées seconde soient à peu près stables. Utilisez `Grid.explicitEuler( D, 2., 0.25)` pour fabriquer un vecteur `E`. 2. il faut calculer le gradient et le Hessien de `E` en tout pixel du domaine. Il suffit de multiplier `E` par l'opérateur gradient selon x `self.Gx` pour obtenir `Ex` et de remultiplier ce vecteur par le même opérateur pour obtenir `Exx`, la dérivée seconde de `E` selon $xx$. Procéder de la même manière pour obtenir `Ey`, `Exy`, `Eyx`, `Eyy`. 3. Pour chaque pixel $(i,j)$, donc pour chaque indice $k$, on va donc réperer si c'est une crête. a. on calcule la matrice Hessienne au pixel d'indice $k$ : $H=\begin{bmatrix} Exx[k] & Exy[k] \\ Eyx[k] & Eyy[k] \end{bmatrix}$ b. on extrait ses valeurs propres $W$ et vecteurs propres $V$, triés dans l'ordre croissant. On a alors dans les notations de [Ridge Detection (Wikipedia)](https://en.wikipedia.org/wiki/Ridge_detection) $\quad L_{pp} = W[0]$, $L_{qq} = W[1]$, qui sont liées aux **courbures principales** c. on fabrique le vecteur $(L_p,L_q)$ comme le gradient $(Ex[k],Ey[k])$ dans la base $(V[:,0],V[:,1])$, $\quad L_p := Ex[k]*V[0,0]+Ey[k]*V[1,0]$ $\quad L_q := Ex[k]*V[0,1]+Ey[k]*V[1,1]$ d. on applique le critère de **ridge** (ou crête) pour notre pixel: $\quad$ **pixel $(i,j)$ est une crête** ssi $|L_p| \le \epsilon$ et $L_{pp} \le v$ et $|L_{pp}| \ge |L_{qq}|$ En pratique, on pourra prendre $\epsilon \approx 0.5$, $v \approx -0.1$. ::: Voilà le genre de résultats que vous pouvez obtenir, ce qui constitue une assez bonne approximation du squelette / axe médian. Pour être propre, il faudrait enlever les quelques "ridges" détectés sur les bords du domaine, car nous avons mal défini le opérateurs dérivatifs à ces endroits. | domaine $\Omega$ | bord de $\Omega$ | Axe médian / squelette | | --------------------------------------------------------------------------------------------------- | -------------------------------------------------------------------------------------------------------- | --------------------------------------------------------------------------------------------------------------- | |  |  |  | |  |  |  | :::success That's all folks ! :::