TI/Programowanie dla Fizyków Medycznych/Dicom

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

"Programowanie dla Fizyków Medycznych"