TI/Programowanie dla Fizyków Medycznych/Obsługa plików graficznych i DICOM

Wczytywanie plików graficznych

Pliki graficzne w większości najpopularniejszych formatów można otworzyć w pythonie korzystając z PIL (Python Image Library). Można to zrobić następującym poleceniem.

from PIL import Image
import numpy as np
import pylab as py
obraz=Image.open('logo2.bmp')
obraz.show()

Jak widać wczytany obiekt nie jest tablicą numpy, ma za to na przykład własną metodę pozwalającą na wyświetlanie wczytanego pliku. Dla naszej wygody możemy jednak przetworzyć tak wczytany plik do postaci tablicy.

im = np.asarray(obraz)
py.imshow(im,cmap = py.cm.gray, interpolation = 'nearest')
py.show()

Obsługa plików DICOM - pyDICOM

Format plików DICOM, czyli Digital Imaging and Communications in Medicine (Obrazowanie Cyfrowe i Wymiana Obrazów w Medycynie) to norma opracowana przez ACR/NEMA (American College of Radiology / National Electrical Manufacturers Association). Do osługi plików tego typu służy między innymi biblioteka pyDICOM. Za jej pomocą wczytujemy przykładowy plik dostępny tutaj Media:I00001.dcm następującym poleceniem

import dicom
plik=dicom.read_file('I00001.dcm')

Plik DICOM poza danymi samego obrazu medycznego zawiera wiele dodatkowych informacji na temat przeprowadzonego badania. Listę dostępnych danych możemy uzyskać poleceniem

print plik.dir()
#>>>['AccessionNumber', 'AcquisitionDate', 'AcquisitionDeviceProcessingCode', 'AcquisitionDeviceProcessingDescription', 'AcquisitionNumber', 'AcquisitionTime', 'AnnotationDisplayFormatID', 'BitsAllocated', 'BitsStored', 'BodyPartExamined', 'BorderDensity', 'Columns', 'ContentDate', 'ContentTime', 'ContrastBolusAgent', 'DerivationDescription', 'FilmOrientation', 'HighBit', 'ImageDisplayFormat', 'ImageType', 'ImagerPixelSpacing', 'InstanceNumber', 'InstitutionName', 'Laterality', 'LossyImageCompression', 'Manufacturer', 'Modality', 'NumberOfPoints', 'PatientBirthDate', 'PatientID', 'PatientName', 'PatientOrientation', 'PatientSex', 'PerformedProcedureStepDescription', 'PerformedProcedureStepID', 'PhotometricInterpretation', 'PixelData', 'PixelRepresentation', 'PixelSpacing', 'PlateID', 'PositionerType', 'RefdStudySequence', 'ReferencedStudySequence', 'ReferringPhysicianName', 'RequestAttributesSequence', 'RequestingService', 'RescaleIntercept', 'RescaleSlope', 'RescaleType', 'Rows', 'SOPClassUID', 'SOPInstanceUID', 'SamplesPerPixel', 'Sensitivity', 'SeriesDate', 'SeriesInstanceUID', 'SeriesNumber', 'SeriesTime', 'SpecificCharacterSet', 'StationName', 'StudyDate', 'StudyDescription', 'StudyID', 'StudyInstanceUID', 'StudyTime', 'Trim', 'ViewPosition']

Podana list zawiera pola które posiada wczytany obiekt. W ten sposób mamy dostęp na przykład do nazwiska pacjenta czy daty badania.

print plik.PatientsName
print plik.AcquisitionDate
#>>>GUBIEC^TOMASZ
#>>>20141119

Dane te możemy dodatkowo edytować.

plik.PatientsName='Jan Kowalaski'
print plik.PatientsName
>>>Jan Kowalaski

Najciekawsze jednak dla nas będą dane opisujące sam obraz. Możemy uzyskać do nich dostęp w wygodnej dla nas postaci tablicy numpy.

pixel=plik.pixel_array
print pixel.shape
print min(pixel.flatten()),max(pixel.flatten())
#py.imshow(pixel,cmap=py.cm.gray,interpolation='nearest')
#py.show()
#>>>(2964, 2364)
#>>> 0 3908

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

"Programowanie dla Fizyków Medycznych"