import numpy as np
import matplotlib.pyplot as plt
def tas_inserer(tas, element):
"""Ajoute element dans le tas"""
tas.append(element) # on place l'élément à la fin
k = len(tas) - 1
# On le fait remonter tant qu'il est plus petit que son parent
while k > 0:
parent = (k - 1) // 2
if tas[k] < tas[parent]:
tas[k], tas[parent] = tas[parent], tas[k] #on les échange
k = parent
else:
break
def tas_extraire_min(tas):
"""Retire et renvoie le plus petit élément du tas"""
racine = tas[0] # on gare la valeur dans la variable racine
dernier = tas.pop() # on retire le dernier élément
n = len(tas)
if n > 0:
tas[0] = dernier # on le place à la racine("écrasement")
k = 0
# On le fait descendre tant qu'un enfant est plus petit
while True:
gauche = 2 * k + 1
droite = 2 * k + 2
petit = k
if gauche < n and tas[gauche] < tas[petit]:
petit = gauche
if droite < n and tas[droite] < tas[petit]:
petit = droite
if petit == k:
break
tas[k], tas[petit] = tas[petit], tas[k]
k = petit
return racine
# ============================================================
# RÉSOLUTION LOCALE (équation du second degré, schéma upwind)
# ============================================================
def temps_local(a, b, f, h):
"""
Calcule T en un point à partir de a et b, les temps des meilleurs
voisins Known
a ou b vaut +inf s'il n'y a pas de voisin Known dans cette direction.
Équation upwind : (T - a)^2 + (T - b)^2 = (h / f)^2
"""
c = h / f # temps pour traverser une case à la vitesse f
# Aucun voisin Known (ne devrait pas arriver)
if a == np.inf and b == np.inf:
return np.inf
# Une seule direction utilisable : équation du premier degré
if a == np.inf or b == np.inf:
return min(a, b) + c
# Les deux directions sont utilisables.
# Équation du second degré : 2 T^2 - 2 (a+b) T + a^2 + b^2 - c^2 = 0
# Le discriminant réduit vaut delta = 2 c^2 - (a - b)^2 > 0 ici
delta = 2 * c**2 - (a - b)**2
return (a + b + np.sqrt(delta)) / 2
def calcul_point(T, etat, F, i, j, h):
"""
Calcule le nouveau temps du point (i, j) Ã partir de ses voisins
Known. Les bords de la grille sont gérés ici.
"""
ny, nx = T.shape
# Meilleur voisin Known dans la direction y (haut et bas)
a = np.inf
for k in [i - 1, i + 1]:
if 0 <= k < ny and etat[k, j] == 2: # on vérifie que le point appartient bien à la grille et que c'est un kn
a = min(a, T[k, j]) # si un des 2 conditions n'est pas respecté, a reste à +inf
# Meilleur voisin Known dans la direction x (gauche et droite)
b = np.inf
for l in [j - 1, j + 1]: # même chose pour b
if 0 <= l < nx and etat[i, l] == 2:
b = min(b, T[i, l])
return temps_local(a, b, F[i, j], h)
# ============================================================
# FAST MARCHING
# ============================================================
def fast_marching(F, sources, h=1.0):
"""
Renvoie la carte T des temps d'arrivée du front.
F : tableau 2D de vitesses (strictement positives)
sources : liste de points de départ [(i, j), ...]
(i = numéro de ligne, j = numéro de colonne)
h : pas de la grille
"""
ny, nx = F.shape # ny lignes, nx colonnes
# ----- Initialisation -----
T = np.ones((ny, nx)) * np.inf # temps +inf = pas encore atteint
etat = np.zeros((ny, nx), dtype=int)
# etat : 0 = Far, 1 = Trial, 2 = Known
tas = [] # tas : éléments (temps, i, j)
voisins = [(-1, 0), (1, 0), (0, -1), (0, 1)]
# Les sources sont Known, de temps nul
for (i0, j0) in sources:
T[i0, j0] = 0.0
etat[i0, j0] = 2
tas_inserer(tas,(0.0,i0,j0))
# ----- Propagation du front -----
while len(tas) > 0:
t, i, j = tas_extraire_min(tas) # on retire le Trial de plus petit temps
# Trial -> Known
etat[i, j] = 2
# Mise à jour des voisins qui ne sont pas encore Known
for (di, dj) in voisins:
ni = i + di
nj = j + dj
if 0 <= ni < ny and 0 <= nj < nx and etat[ni, nj] != 2: # on vérifie que le point est dans la grille et que ce n'est pas un knwon
t_new = calcul_point(T, etat, F, ni, nj, h) # on calcule le nouveau temps
if t_new < T[ni, nj]:
T[ni, nj] = t_new
etat[ni, nj] = 1
tas_inserer(tas, (t_new, ni, nj))
return T
# AFFICHAGE
# ============================================================
def tracer(T, sources, titre, h=1.0, niveaux=12):
"""Carte des temps (Ã gauche) et isochrones (Ã droite)."""
ny, nx = T.shape
x = np.arange(nx) * h
y = np.arange(ny) * h
plt.figure(figsize=(12, 5))
# Carte des temps
plt.subplot(1, 2, 1)
plt.contourf(x, y, T, niveaux)
plt.colorbar(label="Temps T")
for (i0, j0) in sources:
plt.plot(j0 * h, i0 * h, 'r*', markersize=12)
plt.axis('equal')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Carte des temps - " + titre)
# Isochrones
plt.subplot(1, 2, 2)
cs = plt.contour(x, y, T, niveaux)
plt.clabel(cs, fontsize=7)
for (i0, j0) in sources:
plt.plot(j0 * h, i0 * h, 'r*', markersize=12)
plt.axis('equal')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Isochrones - " + titre)
plt.tight_layout()
plt.show()
# ============================================================
# TESTS
# ============================================================
# ---------- Test 1 : vitesse constante F = 1 ----------
ny, nx, h = 100, 100, 1.0
F1 = np.ones((ny, nx))
sources = [(25, 25)]
T1 = fast_marching(F1, sources, h)
tracer(T1, sources, "F = 1", h)
# ----------Test 2: vitesse augmente au fur et à mesure qu'elle s'éloigne de la source -----------
'''ny, nx, h = 100, 100, 1.0
F= np.ones((ny,nx))
sources = [(50,50)]
for i in range(1, ny,1):
for j in range (1, nx,1):
F[i,j] = 10**np.sqrt(((i-ny/2)/100)**2+((j-nx/2)/100)**2)
T = fast_marching(F, sources, h)
tracer(T,sources, "F= 2 fois plus rapide au fur et à mesure qu'on s'éloigne", h)'''
#----------Test 3: vitesse augmente dans vers la gauche mais constante vers la droite ----------
ny, nx, h = 100, 100, 1.0
F= np.ones((ny,nx))
sources = [(50,50)]
for i in range(1, ny,1):
for j in range (1, nx,1):
if j >= 50 :
F[i,j] = 10**np.sqrt(((i-ny/2)/100)**2+((j-nx/2)/100)**2)
T = fast_marching(F, sources, h)
tracer(T, sources, "F=augmente vers vers le haut et constante vers le bas",h )
#-----------Test 4: vitesse aléatoire-----------
ny, nx, h = 100, 100, 1.0
F = 0.5 + np.random.rand(ny, nx) # nombre aléatoire entre 0 et 1
sources = [(50,50)]
T = fast_marching(F, sources, h)
tracer(T, sources, "F=aléatoire",h )
#----------Test 5: Disque lent ----------------
ny, nx, h = 101, 101, 0.01
i_centre, j_centre, rayon = 50, 50, 15 # centre et rayon du disque (en cases)
# Masque : True dans le disque, False ailleurs
masque = np.fromfunction(
lambda i, j: (i - i_centre)**2 + (j - j_centre)**2 < rayon**2, (ny, nx))
F = np.ones((ny, nx))
F[masque] = 0.2 # vitesse faible dans le disque
sources = [(50, 20)] # source à gauche du disque
T = fast_marching(F, sources, h)
tracer(T, sources, "disque lent", h)
#---------Partie 3----------
from PIL import Image
# 1. Charger l'image en niveaux de gris
I = np.array(
Image.open("image.png").convert("L"),
dtype=float
) / 255.0
# 2. Calculer le gradient et sa norme
gy, gx = np.gradient(I)
G = np.sqrt(gx**2 + gy**2)
# 3. Construire la carte de vitesse de la forme 1/(1+a+G)
alpha = 15
F = 1 / (1 + alpha * G)
sources = [(50, 50)]
T = fast_marching(F, sources, h=1.0)
tracer(T, sources, "Segmentation de l'image")
To embed this project on your website, copy the following code and paste it into your website's HTML: