!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')

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()

import os file_size = os.stat('./boats.tif').st_size print(f"La taille de l'image est de", file_size, "octets")

print (img.GetSize()) print (img.GetPixelIDTypeAsString())

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[10:-10,:])

affiche_img(img[10:200,0:200])

affiche_img(img[510:570,340:500])

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 maximale de l'échantillon est de", sub_arr.max()) print(f"La valeur minimale de l'échantillon est de", sub_arr.min()) print(f"L'écart type de l'échantillon est de", sub_arr.std()) print(f"La valeur de la variance de l'échantillon est de", sub_arr.var()) print(f"La valeur moyenne de l'échantillon est de ", sub_arr.mean()) print(f"La somme des valeurs de l'échantillon est de ", sub_arr.sum())

sub_img = sitk.GetImageFromArray(sub_arr)

stats = sitk.StatisticsImageFilter() stats.Execute(sub_img)

print(f"La valeur minimale de l'échantillon est de", stats.GetMinimum()) print(f"La valeur maximale de l'échantillon est de", stats.GetMaximum()) print(f"La valeur moyenne de l'échantillon est de", stats.GetMean()) print(f"La valeur de la variance de l'échantillon est de", stats.GetVariance()) print(f"La valeur de l''ecart type de l'échantillon est de", stats.GetSigma()) print(f"La valeur de la somme de l'échantillon est de", stats.GetSum())

help(np.histogram)

histo = np.histogram(sitk.GetArrayFromImage(img),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)

print (len(histo[0])) print (len(histo[1]))

hist_equal = sitk.AdaptiveHistogramEqualization(img, alpha=10, 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 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ées 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_8bitspoisson = sitk.GetArrayFromImage(img8bit) arr_img8bit = np.random.randint(0,255,4096).reshape((64,64)) arr_8bitspoisson = np.random.poisson(100,4096).reshape((64,64)) hist_img8bit = np.histogram(arr_img8bit) hist_poisson = np.histogram(arr_8bitspoisson) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title('Image générée avec du bruit aléatoire') plt.imshow(arr_img8bit,cmap=plt.cm.Greys_r,interpolation='none') plt.subplot(1,2,2) plt.title('Image générée avec du bruit de Poisson') plt.imshow(arr_8bitspoisson,cmap=plt.cm.Greys_r,interpolation='none')

#On génère l'histogramme de poisson histo_poisson = np.histogram(arr_8bitspoisson,bins =50, density=True) #On génère l'histogramme de poisson cumulé histoCum_poisson = np.cumsum(histo_poisson[0]) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title('Histogramme normalisé poisson') plot = plt.plot(histo_poisson[1][1:],histo_poisson[0]) plt.subplot(1,2,2) plt.title('Histogramme cumulé de poisson') plot = plt.plot(histo_poisson[1][1:],histoCum_poisson)

#On génère l'histogramme de la figure 1 histo_img8bit = np.histogram(arr_img8bit,bins =50, density=True) #On génère l'histogramme de poisson cumulé histoCum_8bit = np.cumsum(histo_img8bit[0]) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title('Histogramme normalisé bruit aléatoire') plot = plt.plot(histo_img8bit[1][1:],histo_img8bit[0]) plt.subplot(1,2,2) plt.title('Histogramme cumulé bruit aléatoire') plot = plt.plot(histo_img8bit[1][1:],histoCum_8bit)

ce_1 = sitk.ReadImage('./CE_1.tif') ce_2 = sitk.ReadImage('./CE_2.tif') arr1 = sitk.GetArrayFromImage(ce_1) arr2 = sitk.GetArrayFromImage(ce_2) plt.subplots(1,2, figsize=(15,12)) plt.subplot(2,2,1) plt.imshow(arr1[0]) plt.subplot(2,2,2) plt.imshow(arr1[0],cmap=plt.cm.Greys_r,interpolation='none') plt.subplot(2,2,3) plt.imshow(arr2[1]) plt.subplot(2,2,4) plt.imshow(arr2[1],cmap=plt.cm.Greys_r,interpolation='none')

#On génère l'histogramme de la figure 1 histo_arr1 = np.histogram(arr1,bins =50, density=True) #On génère l'histogramme cumulé histoCum_arr1 = np.cumsum(histo_arr1[0]) #On génère l'histogramme de la figure 2 histo_arr2 = np.histogram(arr2,bins =50, density=True) #On génère l'histogramme cumulé histoCum_arr2 = np.cumsum(histo_arr2[0]) plt.subplots(1,2,figsize=(15,10)) plt.subplot(2,2,1) plt.title('Histogramme normalisé de l\'image scintigraphique antérieur') plot = plt.plot(histo_arr1[1][1:],histo_arr1[0]) plt.xscale('log') plt.yscale('log') plt.subplot(2,2,2) plt.title('Histogramme cumulé de l\'image scintigraphique antérieur') plot = plt.plot(histo_arr1[1][1:],histoCum_arr1) plt.xscale('log') plt.yscale('log') plt.subplot(2,2,3) plt.title('Histogramme normalisé de l\'image scintigraphique postérieur') plot = plt.plot(histo_arr2[1][1:],histo_arr2[0]) plt.xscale('log') plt.yscale('log') plt.subplot(2,2,4) plt.title('Histogramme cumulé de l\'image scintigraphique postérieur') plot = plt.plot(histo_arr2[1][1:],histoCum_arr2) plt.xscale('log') plt.yscale('log')

help(sitk.HistogramMatching)

im_match = sitk.HistogramMatching(ce_2,ce_1, numberOfMatchPoints=100) im_match_image = sitk.GetArrayFromImage(im_match) plt.imshow(im_match_image[1].T,cmap=plt.cm.Greys) plt.show() histogramatch = np.histogram(sitk.GetArrayFromImage(im_match),bins =50, density=True) histoCummatch = np.cumsum(histogramatch[0]) plt.subplots(1,2,figsize=(15,5)) plt.subplot(1,2,1) plt.title('Histogramme matching normalisé') plt.plot(histogramatch[1][1:],histogramatch[0]) plt.plot(histo_arr1[1][1:],histo_arr1[0]) plt.plot(histo_arr2[1][1:],histo_arr2[0]) plt.legend(['Histogramme matching', 'ce_1', 'ce_2']) plt.yscale('log') plt.subplot(1,2,2) plt.title('Histogramme matching cumulé') plt.plot(histogramatch[1][1:],histoCummatch) plt.plot(histo_arr1[1][1:],histoCum_arr1) plt.plot(histo_arr2[1][1:],histoCum_arr2) plt.legend(['Histogramme matching', 'ce_1', 'ce_2']) plt.yscale('log')

print(f"La valeur minimale de l'échantillon est de", arr1[1].max()) print(f"La valeur maximale de l'échantillon est de", arr2[1].max())

plt.subplots(1,3, figsize=(15,8)) plt.subplot(3,1,1) plt.imshow(arr1[1].T,cmap=plt.cm.Greys,interpolation='none') plt.title('Image postérieur 1') plt.subplot(3,1,2) plt.imshow(arr2[1].T,cmap=plt.cm.Greys,interpolation='none') plt.title('Image postérieur 2') plt.subplot(3,1,3) plt.imshow(im_match_image[1].T,cmap=plt.cm.Greys,interpolation='none') plt.title('Image postérieur postmatching')

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) #plt.axis('off');

img_seuillage = sitk.Threshold(ce_1,0,190) arr_seuillage = sitk.GetArrayFromImage(img_seuillage) plt.figure(figsize=(15,10)) plt.imshow(arr_seuillage[1].T, cmap=plt.cm.Greys) plt.axis('off');

plt.subplots(1,3, figsize=(15,5)) plt.subplot(2,2,1) plt.imshow(sitk.GetArrayFromImage(sitk.OtsuThreshold(ce_1[:,:,1])).T) plt.subplot(2,2,2) plt.imshow(sitk.GetArrayFromImage(sitk.LiThreshold(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) plt.show()

def MaConvolution_kern1(monimage): ker_1 = np.ones((3,3),dtype=np.int8) arr_boat = sitk.GetArrayFromImage(monimage) 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 = (np.ones((3,3),dtype=np.int16)) imageconvoluee = sitk.Convolution(sitk.Cast(img,sitk.sitkInt16),sitk.GetImageFromArray(ker), normalize=True) 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(sitk.GetArrayFromImage(imageconvoluee), cmap=plt.cm.Greys_r,interpolation='none');

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) 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_bateau = 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.imshow(gauss_bateau,cmap=plt.cm.Greys_r,interpolation='none') plt.show()

def MaConvolution_kern3(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

kern_sobel_y = np.zeros((3,3), dtype=np.float32) kern_sobel_x = np.zeros((3,3), dtype=np.float32) kern_sobel_y[0,0] = -1 kern_sobel_y[0,1] = 0 kern_sobel_y[0,2] = 1 kern_sobel_y[1,0] = -2 kern_sobel_y[2,0] = -1 kern_sobel_y[2,2] = 1 kern_sobel_y[1,2] = 2 conv_x = MaConvolution_kern3(img, kern_sobel_y.T) conv_y = MaConvolution_kern3(img, kern_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()

kern_sobel_y = np.zeros((3,3), dtype=np.float32) kern_sobel_y[0,0] = kern_sobel_y[0,2] = -1 kern_sobel_y[0,1] = -2 kern_sobel_y[2,0] = kern_sobel_y[2,2] = 1 kern_sobel_y[2,1] = 2 Sobel_boat_y = sitk.Convolution(sitk.Cast(img,sitk.sitkFloat32), sitk.GetImageFromArray(kern_sobel_y)) Sobel_force = sitk.Sqrt(sitk.Square(Sobel_boat_y) + sitk.Square(Sobel_boat_y)) plt.figure(figsize=(10,10)) plt.imshow(sitk.GetArrayFromImage(Sobel_force),cmap=plt.cm.Greys_r,interpolation='none') plt.show()

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 = k+1 ar_img[h,l] = np.median(mask) plt.figure(figsize=(10,10)) img_median = sitk.GetImageFromArray(ar_img) plt.imshow(ar_img, cmap=plt.cm.Greys_r,interpolation='none') plt.show()

%%time imgDenoise = sitk.Median(imgPoivreSel) plt.subplots(figsize=(15,15)) plt.subplot(1,3,1) plt.imshow(sitk.GetArrayFromImage(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(sitk.GetArrayFromImage(imgDenoise) - ar_img,cmap=plt.cm.Greys_r,interpolation='none')

plt.subplots(figsize=(10,10)) plt.subplot(1,2,1) plt.imshow(sitk.GetArrayFromImage(imgPoivreSel),cmap=plt.cm.Greys_r,interpolation='none') plt.subplot(1,2,2) plt.imshow(ar_img,cmap=plt.cm.Greys_r,interpolation='none') plt.show()