#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Tue Dec 14 18:38:14 2021 @author: maximebouchereau """ # Projet - problème de tomographie import numpy as np import matplotlib.pyplot as plt import random from matplotlib import image import sys # Question 1 def Num(Nx,Ny,j,i): """Donne le numéro du bloc en position (i,j) dans le modèle, la position utilisée étant la même que celle utilisée dans le cas des matrices. Entrées: - (Nx,Ny): Entiers - Nombre de cellulles (Horizontal,Vertical) - (i,j): Entiers - Coordonnées de la cellule sur la grille""" return i + (j-1)*Nx def G(Nx,Ny): """Donne la matrice G à N lignes et K colonnes telle que: - N est le nombre de rayons traversant le modèle avec la géométrie souhaitée dans le projet - K est le nombre de blocs de notre modèle - Pour tous 1 <= j <= N, 1 <= k <= K, l'élémlent d'indice (j,k) de cette matrice correspond à la longueur du rayon j dans le bloc k - La numérotation des rayons est la suivante: > Les rayons horizontaux sont numérotés de 1 à Ny > Les rayons verticaux sont numérotés de Ny+1 à Nx+Ny > Les rayons obliques orientés vers le haut sont numérotés de Nx+Ny+1 à 2(Nx+Ny)-1 > Les rayons obliques orientés vers le bas sont numérotés de 2(Nx+Ny) à 3(Nx+Ny)-2 - La longueur dans chque bloc est la suivante selon le type de rayon: > Pour les deux premières catégories, le rayon est 1 > Pour les deux dernières catégories, le rayon est de sqrt(2) Entrées: - Nx,Ny: entiers - Dimension du domaine (nombre de blocs en x et y)""" # Paramètres calculés K = Nx*Ny # Nombre de blocs N = 3*(Nx+Ny)-2 # Nombre de rayons MatG = np.zeros((0,K)) matG = np.zeros((1,K)) for j in range(N): matG = np.zeros((1,K)) if 1 <= j <= Ny: # Premier groupe de rayons (orienté E) for i in range(1,Nx+1): matG[0,Num(Nx,Ny,j,i)-1] = 1 MatG = np.concatenate((MatG,matG),axis=0) if Ny+1 <= j <= Ny+Nx: # Deuxième groupe de rayons (orienté S) for i in range(1,Ny+1): matG[0,Num(Nx,Ny,i,j-Ny)-1] = 1 MatG = np.concatenate((MatG,matG),axis=0) if Nx+Ny+1 <= j <= 2*(Nx+Ny)-1: # Troisième groupe de rayons (orienté NE) for i in range(0,Nx+Ny): if 1<= j-(Nx+Ny)-i <= Nx and 1<= i+1 <= Ny: matG[0,Num(Nx,Ny,i+1,j-(Nx+Ny)-i)-1] = np.sqrt(2) MatG = np.concatenate((MatG,matG),axis=0) if 2*(Nx+Ny) <= j <= 3*(Nx+Ny)-2: # Quatrième groupe de rayons (orienté SE) for i in range(0,Nx+Ny): if 1 <= i+1 <= Nx and 1 <= j - 2*(Nx+Ny) + 1 - i <= Ny: matG[0,Num(Nx,Ny,j-2*(Nx+Ny)+1-i,Nx-i)-1] = np.sqrt(2) MatG = np.concatenate((MatG,matG),axis=0) return MatG def Pseudo_Inverse(M): """Calcule la matrice pseudo-inverse d'une matrice M via la méthode SVD Entrée: - M: array de taille arbitraire""" return np.linalg.pinv(M, rcond=1e-15, hermitian=False) def Modele(Im): """Transforme l'image Im en un vecteur colonne (vecteur modèle, m dans l'énoncé) Entrée: - Im: array de taille (Nx,Ny) représentant les valeurs attribuées à chaque pixel d'une image""" Nx = Im.shape[1] Ny = Im.shape[0] K = Nx*Ny m = np.zeros(K) for i in range(Ny): for j in range(Nx): m[Num(Nx,Ny,i+1,j+1)-1] = Im[i,j] return m def Image(Nx,Ny,m): """Transforme le vecteur modèle m en tableau 2D prêt à être converti en image Entrées: - Nx,Ny: entiers - Dimensions auxquelles on veut mettre la tableau - m: array - Vecteur de taille K = NxNy""" Im = np.zeros((Ny,Nx)) for i in range(Ny): for j in range(Nx): Im[i,j] = m[Num(Nx,Ny,i+1,j+1)-1] return Im def Affiche(Im,title): """Affiche une image Im sous forme de nuances de gris Entrée: - Im: array contenant des valeurs, les valeurs les plus petites seront affichées en blanc, les plus grandes en noir - title: chaîne de caractères - Titre de l'image""" plt.figure() plt.title(title) plt.imshow(Im,cmap="binary") def Pert_Bruit(M,sigma): """Perturbe les données d'une matrice M par un bruit d'amplitude sigma. Un bruit d'amplitude 1 correspond à une variable gaussienne centrée réduite Entrées: - M: array de taille quelconque - sigma: réel positif - Amplitude du bruit avec lequel les données sont perturbées""" return M + sigma*np.random.normal(0,1,size=M.shape) def Recree_Image(Im,sigma=0): """Fonction qui prend en entrée une image sous forme d'un tableau de données, va recréer des données en envoyant des rayons et en perturbant par un bruit, puis va reconstituer l'image en inversant la matrice G (pseudo-inverse). Affiche cette image Entrées: - Im: array - Tableau de données donnant l'image en noir et blanc (nuances de gris) - sigma: flottant (réel) positif - Amplitude du bruit qui perturbe les données""" Affiche(Im, "Image originale") [Ny,Nx] = Im.shape Im_model = Modele(Im) # Transforme l'image en un vecteur colonne ou chaque composante est associée à la lenteur dans chaque bloc Im_ray = G(Nx,Ny)@Im_model # Envoie les rayons selon les directions données dans l'énoncé Im_ray_noise = Pert_Bruit(Im_ray, sigma) # Perturbation des données mesurées par un bruit Im_calc = Pseudo_Inverse(G(Nx,Ny))@Im_ray_noise # On reconstitue l'image à partir de ce qu'on a mesuré Im_calc_image = Image(Nx,Ny,Im_calc) Affiche(Im_calc_image,"Image recrée par tomographie") # Questions 2 et 3 def ImPic(Nx,Ny): """Génère un pic, i.e. une image blanche dont le seul point noir est le point central Entrées: -Nx,Ny: entiers - Dimensions de l'image""" Pic = np.zeros((Ny,Nx)) Pic[int(Ny/2),int(Nx/2)] = 1 return Pic # Question 4 [Tracé de la lettre M] M = np.zeros((25,25)) M[5:21,5:7] = np.ones((16,2)) M[5:21,18:20] = np.ones((16,2)) for i in range(6,13): M[i,i] = 1 M[i-1,i] = 1 M[i,24-i] = 1 M[i-1,24-i] = 1 # Question 5 def G_Alea(Nx,Ny,p,trace_schema=False): """Reproduit la fonction G donnée en question 1, mais utilise une distribution aléatoire de rayons Entrées: - Nx,Ny: entiers - Dimension du domaine (nombre de blocs en x et y) - p: entier - Nombre de rayons qui seront séléctionnés au hasard parmi tous ceux possibles - trace_schema: True ou False, False par défaut - trace si besoin le schéma des rayons et du modèle Mode de fonctionnement de l'envoi des rayons: Pour l'envoi d'un rayon, selection au hasard d'un angle et d'un point par lequel va passer le rayon, qui est une ligne droite, puis mesure de la longueur du de chaque rayon dans chaque cellule.""" # Paramètres calculés K = Nx*Ny # Nombre de blocs MatG = np.zeros((0,K)) matG = np.zeros((1,K)) print("Construction aléatoire des rayons...") if trace_schema == True: # Bullshit # Trace le modèle (si demandé) plt.figure() plt.title("Schéma rayons-modèle") XX,YY = np.meshgrid(np.arange(0,Nx+1,1), -np.arange(0,Ny+1,1)) plt.scatter(XX,YY) for i in range(p): sys.stdout.write("\r%d "% int(100*(i+1)/p)+"%") sys.stdout.flush() matG = np.zeros((1,K)) theta = -0.5*np.pi + np.pi*np.random.rand() # Angle du rayon selectionné aléatoirement x_ray,y_ray = Nx*np.random.rand(),-Ny*np.random.rand() # On part au hasard a = np.tan(theta) # Equation cartésienne du rayon: y = ax + b, part d'un point b = y_ray - x_ray*np.tan(theta) # (x_ray,y_ray) avec un angle theta # Va parcourir chaque cellule du modèle et va tester si le rayon traverse la cellule, et, si oui # calculer la longueur du rayon dans la cellule en question, et complète la matrice G au besoin if trace_schema == True: #Bullshit # Trace les rayons sur le schema (si demandé) plt.plot(np.array([x_ray,Nx]),np.array([y_ray,a*Nx+b])) plt.text(x_ray,y_ray,str(i)) for ix in range(Nx): for iy in range(Ny): x_cell = ix + 0.5 # Abscisse du point central de la cellule y_cell = -iy - 0.5 # Ordonnée du point central de la cellule Ix = np.array([x_cell-0.5,a*(x_cell-0.5)+b]) # V Point d'intersection entre le rayon et la droite d'équation x = x_cell - 0.5 Jx = np.array([x_cell+0.5,a*(x_cell+0.5)+b]) # W Point d'intersection entre le rayon et la droite d'équation x = x_cell + 0.5 Iy = np.array([(y_cell-0.5-b)/a,y_cell-0.5]) # H Point d'intersection entre le rayon et la droite d'équation y = y_cell - 0.5 Jy = np.array([(y_cell+0.5-b)/a,y_cell+0.5]) # G Point d'intersection entre le rayon et la droite d'équation y = y_cell + 0.5 # Int est la liste des deux points d'intersection éventuels entre le rayon et le bord de la cellule Int = [] # Deux points au plus (sinon zero) sont les points d'intersecrtion entre le rayon et la cellule if y_cell-0.5 <= Ix[1] <= y_cell+0.5: Int = Int + [Ix] if y_cell-0.5 <= Jx[1] <= y_cell+0.5: Int = Int + [Jx] if x_cell-0.5 <= Iy[0] <= x_cell+0.5: Int = Int + [Iy] if x_cell-0.5 <= Jy[0] <= x_cell+0.5: Int = Int + [Jy] if len(Int)==0: # Au cas o Int = 2*[np.zeros(2)] ind = Num(Nx, Ny, iy+1, ix+1) # Numéro de la cellule dont le centre est (x_cell,y_cell) matG[0,ind-1] = np.linalg.norm(Int[1]-Int[0]) MatG = np.concatenate((MatG,matG),axis=0) return MatG def G_Alea_Bis(Nx,Ny,p): """Reproduit la fonction G donnée en question 1, mais utilise une distribution aléatoire de rayons Entrées: - Nx,Ny: entiers - Dimension du domaine (nombre de blocs en x et y) - p: entier - Nombre de rayons qui seront séléctionnés au hasard parmi tous ceux possibles Mode de fonctionnement de l'envoi des rayons: Départ du bord de l'image, puis marche aléatoire en restant sur le pixel ou bien en se déplaçant sur l'un des 8 pixels alentours. Répétition jusqu'à atteinte du bord, constituant l'envoi d'un rayon. Inconvénients: Centre de l'image moins accessible et peu réaliste physiquement.""" # Paramètres calculés K = Nx*Ny # Nombre de blocs MatG = np.zeros((0,K)) matG = np.zeros((1,K)) print("Construction aléatoire des rayons...") for i in range(p): sys.stdout.write("\r%d "% int(100*(i+1)/p)+"%") sys.stdout.flush() matG = np.zeros((1,K)) ix,iy = 1,1 [ix,iy] = random.choice([[np.random.randint(1,Nx+1),random.choice([1,Ny])],[random.choice([1,Nx]),np.random.randint(1,Ny+1)]]) ixx,iyy = ix,iy ind = Num(Nx,Ny,iy,ix) matG[0,ind-1] = 0.5 while 1 <= ixx <= Nx and 1 <= iyy <= Ny: ixx,iyy = ix + np.random.randint(-1,2), iy + np.random.randint(-1,2) ind = Num(Nx,Ny,iy,ix) indd = Num(Nx,Ny,iyy,ixx) #print("Rayon",i,[ix,iy,ind],"->",[ixx,iyy,indd]) if ixx == ix and iyy != iy: matG[0,ind-1] = matG[0,ind-1] + 0.5 if 0 < iyy < Ny+1: #print("Intérieur") matG[0,indd-1] = matG[0,indd-1] + 0.5 if ixx != ix and iyy == iy: matG[0,ind-1] = matG[0,ind-1] + 0.5 if 0 < ixx < Nx+1: #print("Intérieur") matG[0,indd-1] = matG[0,indd-1] + 0.5 if ixx != ix and iyy != iy: matG[0,ind-1] = matG[0,ind-1] + 0.5*np.sqrt(2) if 0< ixx < Nx+1 and 0 < iyy < Ny+1: #print("Intérieur") matG[0,indd-1] = matG[0,indd-1] + 0.5*np.sqrt(2) ix,iy = ixx,iyy MatG = np.concatenate((MatG,matG),axis=0) return MatG def Recree_Image_Alea(Im,p,sigma=0): """Fonction qui prend en entrée une image sous forme d'un tableau de données, va recréer des données en envoyant p rayons aléatoirement et en perturbant par un bruit, puis va reconstituer l'image en inversant la matrice G (pseudo-inverse). Affiche cette image Entrées: - Im: array - Tableau de données donnant l'image en noir et blanc (nuances de gris) - p: entier - Nombre de rayons que l'on envoie pour étudier l'image - sigma: flottant (réel) positif - Amplitude du bruit qui perturbe les données Mode de fonctionnement de l'envoi des rayons: voir fonction G_Alea""" Affiche(Im, "Image originale") [Ny,Nx] = Im.shape Im_model = Modele(Im) # Transforme l'image en un vecteur colonne ou chaque composante est associée à la lenteur dans chaque bloc GG = G_Alea(Nx,Ny,p) Im_ray = GG@Im_model # Envoie les rayons selon les directions données dans l'énoncé Im_ray_noise = Pert_Bruit(Im_ray, sigma) # Perturbation des données mesurées par un bruit Im_calc = Pseudo_Inverse(GG)@Im_ray_noise # On reconstitue l'image à partir de ce qu'on a mesuré Im_calc_image = Image(Nx,Ny,Im_calc) Affiche(Im_calc_image,"Image recrée par tomographie") def Recree_Image_Alea_Bis(Im,p,sigma=0): """Fonction qui prend en entrée une image sous forme d'un tableau de données, va recréer des données en envoyant p rayons aléatoirement et en perturbant par un bruit, puis va reconstituer l'image en inversant la matrice G (pseudo-inverse). Affiche cette image Entrées: - Im: array - Tableau de données donnant l'image en noir et blanc (nuances de gris) - p: entier - Nombre de rayons que l'on envoie pour étudier l'image - sigma: flottant (réel) positif - Amplitude du bruit qui perturbe les données Mode de fonctionnement de l'envoi des rayons: voir fonction G_Alea_Bis""" Affiche(Im, "Image originale") [Ny,Nx] = Im.shape Im_model = Modele(Im) # Transforme l'image en un vecteur colonne ou chaque composante est associée à la lenteur dans chaque bloc GG = G_Alea_Bis(Nx,Ny,p) Im_ray = GG@Im_model # Envoie les rayons selon les directions données dans l'énoncé Im_ray_noise = Pert_Bruit(Im_ray, sigma) # Perturbation des données mesurées par un bruit Im_calc = Pseudo_Inverse(GG)@Im_ray_noise # On reconstitue l'image à partir de ce qu'on a mesuré Im_calc_image = Image(Nx,Ny,Im_calc) Affiche(Im_calc_image,"Image recrée par tomographie") # Test sur une image quelconque name1 = "superman.jpeg" name2 = "smiley.jpeg" name3 = "yoshi.jpeg" name4 = "space_invaders.jpeg" name5 = "space_invader.png" name6 = "pacman.jpg" name7 = "love.jpeg" name8 = "sunglass.jpg" name9 = "laugh.png" name10 = "kim.jpeg" name11 = "paysage.jpeg" name12 = "chat.jpeg" name = name12 def Gray(name): """Convertit une image en couleurs en nuances de gris de format (Ny,Nx) Entrée: - name: Chaîne de caractères - Nom de l'image que l'on veut mettre en nuance de gris (format jpeg ou png)""" img = image.imread(name) # Permet de transformer une image jpeg en tableau de format (Ny,Nx,3), intensité de R,G,B sur chaque pixel gray = np.zeros(img.shape[0:2]) r, g, b = img[:,:,0], img[:,:,1], img[:,:,2] gray = (1/3) * r + (1/3) * g + (1/3) * b gray = 1 - gray/255 # Pour que les valeurs soient comprises entre 0 et 1 et que le 0 correponde au blanc return gray def Flou(Im,n=25): """Réduit le nombtre de pixels d'une image (tableau de données) Entrée: - Im: array - Image de départ - n: Entier: Nombre de pixels sur le côté x de l'image""" [Ny,Nx] = Im.shape Im_bis = Im[::int(Nx/n),::int(Nx/n)] # Diminue l'image en 25x25 return Im_bis def Img(name,n=25): """Créée une image grise à partir d'un fichier jpeg ou png, après l'avoir transformé en nuances de gris puis réduit le nombre de pixels Entrées: - name: Chaîne de caractères - Nom de l'image que l'on veut mettre en nuance de gris (format jpeg ou png) - n: Entier: Nombre de pixels sur le côté x de l'image""" return Flou(Gray(name),n) # Exemple d'application: # Recree_Image_Alea(Img(name11,n=100),p=100)