TI/Programowanie dla Fizyków Medycznych/Obsługa plików graficznych i DICOM
Z Brain-wiki
Wczytwyanie plików graficznych
from PIL import Image
import numpy as np
import pylab as py
obraz=Image.open('logo2.bmp')
obraz.show()
im = np.asarray(obraz)
py.imshow(im,cmap = py.cm.gray, interpolation = 'nearest')
py.show()
Obsługa plików DICOM - pyDICOM
import dicom
import numpy as np
import pylab as py
plik=dicom.read_file('2.dcm')
pixel=plik.pixel_array
print min(pixel.flatten()),max(pixel.flatten())
pix=np.zeros(pixel.shape,)
pix=np.round(255.0*pixel/max(pixel.flatten()))
print min(pix[:100,:100].flatten()),max(pix[:100,:100].flatten())
kolor=np.zeros((pix.shape[0],pix.shape[1],3))
kolor[:,:,0]=pix
kolor[:,:,1]=pix
kolor[:,:,2]=pix
kolor[:300,:300,0]=1
py.imshow(kolor,interpolation='nearest')
py.show()
Zadanie 4
# -*- coding: utf-8 -*-
# Programowanie dla Fizyków Medycznych
# Mateusz Sitarz
# Plik należy otworzyć w folderze z rozpakowanymi plikami instalacyjnymi DICOM.
import dicom
import numpy
import pylab
##################
# zadanie domowe #
##################
# (zadanie 4 z 2 serii zadań domowych)
'''
Ściągnij plik:
www.fuw.edu.pl/~tgubiec/CDRAD_sample.dcm
Plik wynikowy zapisz pod nazwą będącą pierwotną nazwą badanego pacjenta.
1) Zmień nazwę pacjenta na swoje imię i nazwisko, a datę badania na
swoje urodziny w 2012 roku,
2) usuń szum filtrem medianowym,
3) filtrem gamma popraw widoczność kołek wewnątrz kratki,
4) wybiel tło kółek (kratka może zniknąć) poprzez zamianę wszystkich
wartości powyżej pewnego progu na 255. Wartość progu łatwiej znaleźć
analizując histogram obrazka (pylab.hist(macierz_z_obrazkiem)).
5) Wypełnij małe kółka poprzez zamknięcie,
5) podaj średnicę dużego centralnego koła w centymetrach.
'''
import dicom # TEGO BRAKOWAŁO
# wczytanie pliku:
plan = dicom.read_file("CDRAD_sample.dcm")
print plan.PatientsName
# zamiana nazwy pacjenta:
pierwotna_nazwa = plan.PatientsName
nowa_nazwa = "MateuszSitarz"
plan.PatientsName = nowa_nazwa
# zamiana daty badania:
plan.StudyTime = "20121008"
# wczytanie obrazu DICOM:
obraz = plan.pixel_array
obraz=obraz[800:1200,900:1300]
#pylab.subplot(2,3,1)
#pylab.title(u"obraz")
pylab.imshow(obraz, cmap = pylab.cm.gray, interpolation = 'nearest')
pylab.show()
# ==============================================================================
# filtr medianowy:
def mediana(obraz):
kopia = obraz.copy()
wymiar = kopia.shape
kopia = kopia * 1.0
for x in range(1, wymiar[0] - 1): # (filtru medianowego nie stosuje się dla brzegowych pikseli obrazu)
for y in range(1, wymiar[1] - 1):
otoczenie_piksela = kopia[x-1:x+2, y-1:y+2] # pobranie otoczenia piksela
kopia[x,y] = numpy.median(otoczenie_piksela) #TUTAJ RACZEJ MEDIAN zamiast mean
return kopia
obraz = mediana(obraz)
#pylab.subplot(2,3,2)
#pylab.title(u"filtr medianowy")
#pylab.imshow(obraz_medianowy, cmap = pylab.cm.gray, interpolation = 'nearest')
# ==============================================================================
# filtr gamma:
def gamma(obraz, g = 0.4):
kopia = obraz.copy()
wymiar = kopia.shape
kopia = kopia * 1.0
maximum = numpy.max(kopia)
for x in range(0, wymiar[0]):
for y in range(0, wymiar[1]):
kopia[x,y] = kopia[x,y] / maximum # (normalizacja obrazu)
kopia[x,y] = kopia[x,y] ** g # (zastosowanie definicji operacji gamma - potrzebna normalizacja)
kopia[x,y] = kopia[x,y] * maximum # (powrót do normalnych wartości obrazu)
return kopia
#obraz_gamma = gamma(obraz)
#pylab.subplot(2,3,3)
#pylab.title(u"filtr gamma")
#pylab.imshow(obraz_gamma, cmap = pylab.cm.gray, interpolation = 'nearest')
# ==============================================================================
#pylab.hist(obraz)
#pylab.show()
# wybielanie tła:
#obraz_bialy = obraz.copy()
prog=13800 #wybrane na podstawie histogramu
#fragment_obrazu = obraz_bialy[950:1250, 60:120]
maxim = numpy.max(obraz)
for x in range(0, obraz.shape[0]):
for y in range(0, obraz.shape[1]):
if obraz[x,y]<prog: obraz[x,y]=0
else: obraz[x,y]=255
#for n, val in enumerate(obraz_bialy.flat):
# if val >= prog: # (koła są jaśniejsze od tła)
# obraz_bialy.flat[n] = 255
#pylab.subplot(2,3,4)
#pylab.title(u"wybielone tlo")
pylab.imshow(obraz, cmap = pylab.cm.gray, interpolation = 'nearest')
pylab.show()
# ...
# ==============================================================================
# operacja zamknięcia (erozja, następnie dylacja):
def erosion(obrazek, pedzel):
wynik = obrazek.copy()
wym_pedzla = numpy.shape(pedzel)
pol_wym_pedzla = wym_pedzla[0] / 2
for i in range(0 + pol_wym_pedzla, obrazek.shape[0] - pol_wym_pedzla):
for j in range(0 + pol_wym_pedzla, obrazek.shape[1] - pol_wym_pedzla):
lista = []
for k in range(0, pedzel.shape[0]):
for l in range(0, pedzel.shape[1]):
if pedzel[k,l]:
lista.append(obrazek[i - pedzel.shape[0] + k, j - pedzel.shape[1] + l])
wynik[i,j] = numpy.min(lista)
return wynik
def dilation(obrazek, pedzel):
wynik = obrazek.copy()
wym_pedzla = numpy.shape(pedzel)
pol_wym_pedzla = wym_pedzla[0] / 2
for i in range(0 + pol_wym_pedzla, obrazek.shape[0] - pol_wym_pedzla):
for j in range(0 + pol_wym_pedzla, obrazek.shape[1] - pol_wym_pedzla):
lista = []
for k in range(0, pedzel.shape[0]):
for l in range(0, pedzel.shape[1]):
if pedzel[k,l]:
lista.append(obrazek[i - pedzel.shape[0] + k, j - pedzel.shape[1] + l])
wynik[i,j] = numpy.max(lista)
return wynik
def rysuj_pedzel(a):
'''Rysowanie kolowego pedzla w macierzy kwadratowej (2*a+1) * (2*a+1) (nieparzyste wymiary - SRODEK W 1 PIKSELU). Pedzel wewnatrz kolka ma wartosc 1, a poza nim 0.'''
if a == 1:
tlo = numpy.array([[0,1,0],[1,1,1],[0,1,0]])
if a > 1:
n = 2 * a + 1
tlo = numpy.zeros([n, n])
r = n/2.0
for x in range(-n/2, n/2 + 1):
for y in range(-n/2, n/2 + 1):
if ( (x**2 + y**2) <= r**2 ) :
tlo[x + n/2][y + n/2] = 1
return tlo
moj_pedzel = rysuj_pedzel(1)
obraz_erozja = erosion(obraz, moj_pedzel)
obraz_dylacja_erozja = dilation(obraz_erozja, moj_pedzel)
#pylab.subplot(2,3,5)
#pylab.title(u"zamkniecie")
pylab.imshow(obraz_dylacja_erozja, cmap = pylab.cm.gray, interpolation = 'nearest')
pylab.show()
# ==============================================================================
# średnica centralnego koła:
skala = plan.PixelSpacing
#TUTAJ NIE BARDZO ROZUMIEM CO SIE DZIEJE - ZA PEWNE POTRZEBNA JEST JAKAS NIEWIELKA ZMIANA
# stworzenie obrazu czarno-białego:
cien = obraz.copy()
fragment_cienia = cien[950:1250, 60:120]
szum_cienia = numpy.max(fragment_cienia)
prog_cienia = 1.0*szum_cienia
for n, val in enumerate(cien.flat):
if val < prog_cienia:
cien.flat[n] = 0
elif val >= prog_cienia:
cien.flat[n] = 1
#pylab.subplot(2,3,6)
#pylab.title(u"obraz czarno-bialy")
#pylab.imshow(cien, cmap = pylab.cm.gray, interpolation = 'nearest')
# ...
# szukanie średnicy na OY:
# (kiedy już tło jest białe, a koła czarne)
lista_y = []
for y in range(900, 1300): # skrócenie granic, aby obliczać tylko dla obszaru z największym kołem
liczba_jedynek = 0
for x in range(800, 1200): # skrócenie granic, aby obliczać tylko dla obszaru z największym kołem
if cien[x, y] == 1:
liczba_jedynek = liczba_jedynek + 1
lista_y.append(liczba_jedynek)
max_y = numpy.max(lista_y) # najszersza kolumna - jest średnicą
# zamiana pikseli na milimetry:
y_mm = max_y * skala[1]
print "srednica najwiekszego kola:", y_mm, "mm"
# ==============================================================================
#pylab.show()
# zapisanie pliku (jako pierwotna nazwa pacjenta):
#plan.save_as("%(0).4f.dcm" %{"0":pierwotna_nazwa})
