!pip install SimpleITK==2.1.1

%matplotlib inline import SimpleITK as sitk import numpy as np import matplotlib.pyplot as plt

img = sitk.ReadImage('./boats.tif')

print (img.GetSize()) print (img.GetOrigin()) print (img.GetSpacing()) print (img.GetDirection()) print (img.GetNumberOfComponentsPerPixel())

def affiche_img(image): plt.figure(figsize=(20,10)) arr_img = sitk.GetArrayFromImage(image) plt.imshow(arr_img, cmap=plt.cm.Greys_r,interpolation='none'); plt.axis('on'); plt.show()

affiche_img(img) # plt.figure(figsize=(20,10)) # arr_img = sitk.GetArrayFromImage(img) # plt.imshow(arr_img, # cmap=plt.cm.Greys_r,interpolation='none'); # plt.axis('on'); # plt.show()

print (img.GetSize()) print (img.GetNumberOfPixels())

affiche_img(img[10:-10, :]) # on supprime les bandes noires ci-dessus en commencant à 10 au lieu # de 0 et on enlève 10 par "-10" en pas de 1 (par défaut sauf # mention contraire)

affiche_img(img[515:555, 360:500]) #on zoom sur une zone de l'image en définissant les intervalles # selon les lignes (x) et les colonnes (y)

affiche_img(img[15:180, 2:150])

arr_img = sitk.GetArrayFromImage(img) sub_arr = arr_img[2:80,15:200] plt.imshow(sub_arr, cmap=plt.cm.Greys_r);

print(f"La valeur minimale de l'échantillon est ",sub_arr.min()) print(f"La valeur maximale de l'échantillon est ",sub_arr.max()) print(f"La valeur médiane de l'échantillon est ",len(sub_arr)/2) print(f"La valeur moyenne de l'échantillon est ",sub_arr.mean()) print(f"La variance de l'échantillon est ",sub_arr.var()) print(f"L'écart-type de l'échantillon est ",sub_arr.std()) print(f"La somme des valeurs de l'échantillon est ",sub_arr.sum())

sub_img = sitk.GetImageFromArray(sub_arr)

stats = sitk.StatisticsImageFilter() stats.Execute(sub_img) max(sub_img), min(sub_img)

help(np.histogram)

histo = np.histogram(sitk.GetArrayFromImage(img),bins=50, density=True) len(histo), len(histo[0]), len(histo[1]) # pour construire 50 histogrammes il faut 51 valeurs # c'est pour cette raison le "51"

histo = np.histogram(sitk.GetArrayFromImage(img),bins = 50, density=True) #si on écit "False" #alors ordonnée fausse histoCum = np.cumsum(histo[0]) #non histo[1] car x et y doivent être de même dimension plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title("Histogramme normalisé") plot = plt.plot(histo[1][1:],histo[0]) plt.subplot(1,2,2) plt.title("Histogramme cumulé") plot = plt.plot(histo[1][1:],histoCum)

# hist_equal = sitk.AdaptiveHistogramEqualization(img, alpha=0, beta=1) hist_equal = sitk.AdaptiveHistogramEqualization(img, alpha=0, beta=1) # On recupère le tableau au format numpy np_histEqual = sitk.GetArrayFromImage(hist_equal) # AFFICHAGE de l'image égalisée plt.subplots(1,2, figsize=(15,5)) plt.subplot(1,2,1) plt.imshow(arr_img, cmap=plt.cm.Greys_r,interpolation='none') plt.title('original') plt.axis('off') plt.subplot(1,2,2) plt.imshow(np_histEqual, cmap=plt.cm.Greys_r,interpolation='none') plt.title("egalisation d'histogramme") plt.axis('off') plt.show()

# On génère l'histogramme histo2 = np.histogram(np_histEqual,bins = 50,density=True) # On génère l'histogramme cumulé histoCum2 = np.cumsum(histo2[0]) plt.subplots(1,2, figsize=(15,5)) plt.subplot(1,2,1) plt.title(u'Histogrammes normalisés avant \n et après égalisation'); plt.plot(histo2[1][1:],histo[0],histo2[1][1:],histo2[0]); #plt.ylabel(u'distribution normalisée'); plt.xlabel(u'intensité des pixels'); plt.ylabel(u'distribution normalisée'); plt.legend((u'avant',u'après'),loc = 2) plt.subplot(1,2,2) plt.title(u'Histogrammes cumulés avant \n et après égalisation'); plt.plot(histo[1][1:],histoCum, histo2[1][1:],histoCum2,'r-+'); plt.xlabel(u'intensité des pixels'); plt.legend((u'avant',u'après'),loc=2) plt.subplots_adjust(wspace = 0.5) plt.show();

img8bit = sitk.Image((64,64), sitk.sitkUInt8) arr_img8bit = sitk.GetArrayFromImage(img8bit) arr_img8bit = np.random.randint(0,255,4096).reshape((64,64)) #arr_poisson = np.random.poisson(100,4096).reshape((64,64)) hist_img8bit = np.histogram(arr_img8bit) #hist_poisson = np.histogram(arr_poisson) plt.imshow(arr_img8bit, cmap=plt.cm.Greys_r,interpolation='none') # QR code #plt.axis('off') histo = np.histogram(sitk.GetArrayFromImage(img8bit),bins = 50, density=True) histoCum = np.cumsum(histo[0]) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title("Histogramme normalisé") plot = plt.plot(histo[1][1:],histo[0]) plt.subplot(1,2,2) plt.title("Histogramme cumulé") plot = plt.plot(histo[1][1:],histoCum)

ce_1 = sitk.ReadImage('./CE_1.tif') ce_2 = sitk.ReadImage('./CE_2.tif')

#print(ce_1.GetSize()) #print(ce_2.GetSize()) ar_ce_1 = sitk.GetArrayFromImage(ce_1) ar_ce_2 = sitk.GetArrayFromImage(ce_2) #print(sitk.GetArrayFromImage(ce_1).shape) plt.figure(figsize=(30,10)) plt.imshow(ar_ce_1[0].T,cmap=plt.cm.Greys) #image antérieure plt.show() plt.figure(figsize=(30,10)) plt.imshow(ar_ce_1[1].T,cmap=plt.cm.Greys) #image postérieure plt.show()

histo1 = np.histogram(ar_ce_1[1],bins = 50, density=True) histoCum1 = np.cumsum(histo1[0]) histo2 = np.histogram(ar_ce_2[1],bins = 50, density=True) histoCum2 = np.cumsum(histo2[0]) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,1,1) plt.title("Histogramme normalisé") plt.plot(histo1[1][1:],histo1[0], label="ce_1") plt.plot(histo2[1][1:],histo2[0], label="ce_2") plt.legend() plt.xscale("log") plt.yscale("log") plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,1,1) plt.title("Histogramme cumulé") plt.plot(histo1[1][1:],histoCum1, label="ce_1") plt.plot(histo2[1][1:],histoCum2, label="ce_2") plt.legend() plt.xscale("log") plt.yscale("log")

help(sitk.HistogramMatching)

im_match = sitk.HistogramMatching(ce_2,ce_1) arr_im_match = sitk.GetArrayFromImage(im_match) #print(sitk.GetArrayFromImage(im_match).shape) plt.figure(figsize=(30,10)) plt.subplot(1,3,1) plt.imshow(arr_im_match[1].T,cmap=plt.cm.Greys) plt.figure(figsize=(30,10)) plt.subplot(1,3,2) plt.imshow(ar_ce_1[1].T,cmap=plt.cm.Greys) plt.figure(figsize=(30,10)) plt.subplot(1,3,3) plt.imshow(ar_ce_2[1].T,cmap=plt.cm.Greys) plt.show()

histo_match = np.histogram(arr_im_match, bins = 50, density=True) histoCum_match = np.cumsum(histo_match[0]) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title("Histogramme normalisé") plot = plt.plot(histo1[1][1:],histo1[0], label="ce_1") plt.plot(histo2[1][1:], histo2[0], label="ce_2") plt.plot(histo_match[1][1:], histo_match[0], label="matching") plt.legend() plt.xscale("log") plt.yscale("log") plt.subplot(1,2,2) plt.title("Histogramme cumulé") plt.plot(histo1[1][1:],histoCum1, label='ce_1') plt.plot(histo2[1][1:],histoCum2, label="ce_2") plt.plot(histo_match[1][1:],histoCum_match, label="matching") plt.xscale("log") plt.yscale("log") plt.legend() #ar_ce_1[1].max(),ar_ce_2[1].max()

arr_post_1 = sitk.GetArrayFromImage(ce_1[:,:,1]) arr_post_1[arr_post_1>190]=0 plt.figure(figsize=(15,10)) plt.imshow(arr_post_1.T, cmap=plt.cm.Greys_r) plt.axis('off');

plt.figure(figsize=(30,10)) plt.imshow(sitk.GetArrayFromImage(sitk.OtsuThreshold(ce_1[:,:,1])).T, cmap=plt.cm.Greys_r)

li_filter = sitk.LiThresholdImageFilter() li_filter.SetInsideValue(0) li_filter.SetOutsideValue(1) li_img = li_filter.Execute(ce_1[:,:,1]) plt.imshow(sitk.GetArrayFromImage(sitk.Cast(li_img,sitk.sitkUInt16)*ce_1[:,:,1]).T, cmap=plt.cm.Greys) sitk.GetArrayFromImage(li_img).T

im_seuillage = sitk.Threshold(ce_1,0,190) arr_seuillage = sitk.GetArrayFromImage(im_seuillage) plt.figure(figsize=(15,10)) plt.imshow(arr_seuillage[1].T, cmap=plt.cm.Greys_r) plt.show()

def MaConvolution_kern1(monimage): #ker_1 = np.ones((3,3),dtype=np.int16) ker_1 = np.ones((3,3),dtype=np.int8) arr_boat = sitk.GetArrayFromImage(monimage) #out_arr = np.array(arr_boat,dtype=np.int16) out_arr = np.array(arr_boat,dtype=np.int8) for u in range(1,arr_boat.shape[0]-1): for v in range(1,arr_boat.shape[1]-1): somme = 0 for i in range(-1,2): for j in range(-1,2): somme = somme + arr_boat[u+i,v+j] out_arr[u,v] = somme return out_arr

%time smooth_img = MaConvolution_kern1(img) plt.subplots(1,2,figsize=(15,15)) plt.subplot(1,2,1) plt.imshow(sitk.GetArrayFromImage(img), cmap=plt.cm.Greys_r,interpolation='none'); plt.subplot(1,2,2) plt.imshow(smooth_img, cmap=plt.cm.Greys_r,interpolation='none');

help(sitk.Convolution)

ker = sitk.GetImageFromArray(np.ones((3,3),dtype=np.int16)) conv_image = sitk.Convolution(sitk.Cast(img,sitk.sitkInt16),ker) #img formatée à 16 bits signés arr_conv_image = sitk.GetArrayFromImage(conv_image) plt.subplots(1,2,figsize=(15,10)) plt.subplot(1,2,1) plt.imshow(sitk.GetArrayFromImage(img),cmap=plt.cm.Greys_r,interpolation='none'); plt.subplot(1,2,2) plt.imshow(arr_conv_image, cmap=plt.cm.Greys_r,interpolation='none');

#création filtre gaussien def gaussienne(x,y,var): return np.exp(-(x*x+y*y)/(2*var))/(2*np.pi*var) dim = 11 kern_Gauss = np.zeros((dim,dim), np.float32) #kernel à 32 bits réels donc img au même format for i in range(-dim//2,dim//2+1): for j in range(-dim//2,dim//2+1): kern_Gauss[i+dim//2,j+dim//2]= gaussienne(i,j,1) gauss_boat = sitk.GetArrayFromImage(sitk.Convolution(sitk.Cast(img,sitk.sitkFloat32), sitk.GetImageFromArray(kern_Gauss),normalize=True)) plt.subplots(1,2,figsize=(15,15)) plt.subplot(1,2,1) plt.title('noyau gaussien') plt.imshow(kern_Gauss,interpolation='none') plt.subplot(1,2,2) plt.title('filtrage gaussien') plt.imshow(gauss_boat,cmap=plt.cm.Greys_r,interpolation='none') plt.show()

help(sitk.DiscreteGaussian)

gauss_im = sitk.DiscreteGaussian(img,2) plt.subplots(1,2,figsize=(10,10)) plt.subplot(1,2,1) plt.title('Filtrage Gaussien par SimpleITK') plt.imshow(sitk.GetArrayFromImage(gauss_im), cmap=plt.cm.Greys_r,interpolation='none') plt.subplot(1,2,2) plt.title('filtrage gaussien analytique') plt.imshow(gauss_boat,cmap=plt.cm.Greys_r,interpolation='none') plt.show()

def MaConvolution_kernel(monimage, kernel): if kernel.shape[0]%2 ==0: raise TypeError ('Le noyau doit être de dimension impaire') else: kx = int(kernel.shape[0]/2) ky = int(kernel.shape[1]/2) ker_som = kernel.sum() arr_boat = sitk.GetArrayFromImage(monimage) out_arr = np.zeros((arr_boat.shape),dtype=np.float32) for v in range (kx, arr_boat.shape[1]-kx): for u in range (ky, arr_boat.shape[0]-ky): somme = 0. for i in range (-kx, kx+1): for j in range (-ky,ky+1): somme += arr_boat[u+i,v+j]*kernel[i+kx,j+ky] if ker_som == 0: out_arr[u,v] = somme else: out_arr[u,v] = somme/ker_som return out_arr

#initialisation des matrice en les posant égales à l'idendité kernel_sobel_y = np.zeros((3,3),dtype=np.float32) kernel_sobel_x = np.zeros((3,3),dtype=np.float32) #création de la matrice Kx de Mr.Sobel kernel_sobel_y[0,0] = -1 kernel_sobel_y[0,1] = 0 kernel_sobel_y[0,2] = 1 kernel_sobel_y[1,0] = -2 kernel_sobel_y[2,0] = -1 kernel_sobel_y[2,2] = 1 kernel_sobel_y[1,2] =2 # les matrices kx et ky sont les transposées de l'une de l'autre conv_x = MaConvolution_kernel(img, kernel_sobel_y.T) conv_y = MaConvolution_kernel(img, kernel_sobel_y)

fig, ax = plt.subplots(1,2,figsize=(15,15)) ax[0].imshow(conv_x, cmap=plt.cm.Greys_r,interpolation='none') ax[1].imshow(conv_y, cmap=plt.cm.Greys_r,interpolation='none') ax[0].axis('off') ax[1].axis('off') plt.show()

#Calcul du module du gradient de l'image sobel_boat_y =sitk.Convolution(sitk.Cast(img,sitk.sitkFloat32), sitk.GetImageFromArray(kernel_sobel_y)) sobel_boat_x =sitk.Convolution(sitk.Cast(img,sitk.sitkFloat32), sitk.GetImageFromArray(kernel_sobel_y.T)) sobel_force =sitk.Sqrt(sitk.Square(sobel_boat_x)+sitk.Square(sobel_boat_y)) arr_sobel_force = sitk.GetArrayFromImage(sobel_force) plt.figure(figsize=(10,10)) plt.imshow(arr_sobel_force, cmap=plt.cm.Greys_r, interpolation='none') plt.show() sobel_edge = sitk.SobelEdgeDetection(sitk.Cast(img,sitk.sitkFloat32) #arr_sob = sitk.GetArrayFromImage(sobel_edge) plt.figure(figsize=(10,10)) plt.imshow(sitk.GetArrayFromImage(sobel_edge), cmap=plt.cm.Greys_r, interpolation='none') plt.show() plt.figure(figsize=(10,10)) plt.imshow(arr_sobel_force - sitk.GetArrayFromImage(sobel_edge), cmap=plt.cm.Greys_r, interpolation='none') plt.show()

help(sitk.SobelEdgeDetection)

imgPoivreSel = sitk.SaltAndPepperNoise(img, 0.05) plt.figure(figsize=(10,10)) plt.imshow(sitk.GetArrayFromImage(imgPoivreSel), cmap=plt.cm.Greys_r,interpolation='none') plt.show()

largeur = imgPoivreSel.GetWidth() hauteur = imgPoivreSel.GetHeight() ar_img = sitk.GetArrayFromImage(imgPoivreSel) mask = np.zeros(9) for l in range(1,largeur-1): for h in range(1,hauteur-1): k = 0 for i in range (-1,2): for j in range (-1,2): mask[k]=ar_img[h+i,l+j] k += 1 ar_img[h,l] = np.median(mask) plt.figure(figsize=(10,10)) plt.imshow(ar_img, cmap=plt.cm.Greys_r,interpolation='none') img_median = sitk.GetImageFromArray(ar_img)

#version plus rapide de l'algorithme précédent imgDenoise = sitk.Median(imgPoivreSel) arr_imgDenoise = sitk.GetArrayFromImage(imgDenoise) plt.subplots(figsize=(10,10)) plt.subplot(1,3,1) plt.imshow(arr_imgDenoise,cmap=plt.cm.Greys_r, interpolation='none') plt.subplot(1,3,2) plt.imshow(ar_img, cmap=plt.cm.Greys_r,interpolation='none') plt.subplot(1,3,3) plt.imshow(arr_imgDenoise - ar_img, cmap=plt.cm.Greys_r, interpolation='none');