from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt


#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
#Etape 1: Chargement et visualisation
# 1. Charger l'image GeoTIFF
filename = "fig_bis/00492.tif"
dataset = gdal.Open(filename)

# 2. Extraire la première bande
band1 = dataset.GetRasterBand(1)
image_array = band1.ReadAsArray()

# 3. Vérifier les dimensions
print("Dimensions de l'image :", image_array.shape)

# Afficher l'image
plt.imshow(image_array, cmap='gray')
plt.title("Une coupe de l'échantillon cylindrique")
plt.colorbar()
plt.show()

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
#Etape 2 : Chargement et visualisation;
# 1 — Segmenter le fond
seuil_fond = np.percentile(image_array, 5)
masque_fond = image_array <= seuil_fond

# 2 — Inverser le masque pour obtenir la roche
masque_roche = 1 - masque_fond.astype(int)
#####
#plt.imshow(masque_fond, cmap="gray")
#plt.show()
#####

# 3 — Trouver les pixels de la roche
coords_roche = np.column_stack(np.where(masque_roche))

# 4 — Calculer le centre
centre = coords_roche.mean(axis=0)

# 5 — Calculer le rayon
distances_roche = np.sqrt((coords_roche[:, 0] - centre[0])**2 + (coords_roche[:, 1] - centre[1])**2)
rayon = distances_roche.max()

print(f"Centre détecté : [{centre[0]:.4f}, {centre[1]:.4f}]")
print(f"Rayon détecté : {rayon:.4f}")

# Affichage final du cylindre détecté
plt.figure(figsize=(6, 6))
plt.imshow(masque_roche, cmap="gray")
#plt.imshow(image_array, cmap="gray")
plt.plot(centre[1], centre[0], "ro")
circle = plt.Circle((centre[1], centre[0]), rayon, color='r', fill=False, linewidth=2)
plt.gca().add_patch(circle)
plt.title("Centre et rayon du cylindre")
plt.axis("off")
plt.show()
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# Étape 3 — Détection des fissures
from scipy.ndimage import uniform_filter, binary_opening
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

seuil_fissure = np.percentile(image_array, 60)
masque_fissure = image_array > seuil_fissure
#masque_fissure = 1 - masque_fissure.astype(int)

plt.figure(figsize=(8, 8))
plt.imshow(masque_fissure, cmap='Blues')
plt.title("Fissures détectées sans marge au bord")
plt.axis("off")
plt.show()
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# 1 — Masque intérieur
marge_bord = 20
y_indices, x_indices = np.indices(image_array.shape)
distances_image = np.sqrt((y_indices - centre[0])**2 + (x_indices - centre[1])**2)
masque_interieur = distances_image < (rayon - marge_bord)

# 2 — Seuil fissure
pixels_roche = image_array[masque_interieur]
seuil_fissure = np.percentile(pixels_roche, 20)
masque_fissure = (image_array < seuil_fissure) & masque_interieur

# 3 — Affichage final
plt.figure(figsize=(8, 8))
plt.imshow(masque_fissure, cmap='Blues')
#plt.imshow(image_rgba, interpolation="none")
plt.title("Fissures détectées avec marge au bord")
plt.axis("off")
plt.show()
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# Etape Optionelle
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# — Nettoyage simple : flou + ouverture morphologique
masque_fissure_flou = uniform_filter(masque_fissure.astype(float), size=3)
masque_fissure = masque_fissure_flou > 0.5

# Ouverture morphologique pour éliminer les petits pixels isolés**
masque_fissure = binary_opening(masque_fissure, structure=np.ones((3, 3)))
hauteur, largeur = image_array.shape
#  Image RGBA fissures bleues
image_rgba = np.zeros((hauteur, largeur, 4), dtype=np.uint8)
image_rgba[masque_fissure] = [0, 0, 255, 255]

# 5 — Affichage final
plt.figure(figsize=(8, 8))
plt.imshow(masque_fissure, cmap='Blues')
#plt.imshow(image_rgba, interpolation="none")
plt.title("Fissures détectées")
plt.axis("off")
plt.show()