TI/Programowanie dla Fizyków Medycznych/Obsługa plików graficznych i DICOM: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 1 pośredniej wersji utworzonej przez tego samego użytkownika)
Linia 16: Linia 16:
  
 
==Obsługa plików DICOM - pyDICOM==
 
==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
+
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 obsł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
 
<source lang="python">
 
<source lang="python">
 
import dicom
 
import dicom
Linia 25: Linia 25:
 
<source lang="python">
 
<source lang="python">
 
print plik.dir()
 
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']
+
>>> ['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']
 
</source>
 
</source>
−
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.
+
Podana lista zawiera pola, które posiada wczytany obiekt. W ten sposób mamy dostęp na przykład do nazwiska pacjenta czy daty badania.
 
<source lang="python">
 
<source lang="python">
 
print plik.PatientsName
 
print plik.PatientsName
Linia 45: Linia 45:
 
print pixel.shape
 
print pixel.shape
 
print min(pixel.flatten()),max(pixel.flatten())
 
print min(pixel.flatten()),max(pixel.flatten())
−
#py.imshow(pixel,cmap=py.cm.gray,interpolation='nearest')
+
py.imshow(pixel,cmap=py.cm.gray,interpolation='nearest')
−
#py.show()
+
py.show()
−
#>>>(2964, 2364)
+
>>>(2964, 2364)
−
#>>> 0 3908
+
>>> 0 3908
 
</source>
 
</source>
 
[[Plik:dicom1.png]]
 
[[Plik:dicom1.png]]
−
 
−
 
−
 
−
 
−
 
−
<source lang="python">
 
−
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()
 
−
</source>
 
−
Zadanie 4
 
−
<source lang="python">
 
−
 
−
# -*- 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})
 
−
</source>
 
  
 
[["Programowanie dla Fizyków Medycznych"]]
 
[["Programowanie dla Fizyków Medycznych"]]

Aktualna wersja na dzień 12:16, 9 cze 2015

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 obsł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 lista 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

Dicom1.png

"Programowanie dla Fizyków Medycznych"

Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=TI/Programowanie_dla_Fizyków_Medycznych/Obsługa_plików_graficznych_i_DICOM&oldid=3612”