« Jarak Titik ke Garis pada bidang 2D
(Another) attempt to volume (3d scalar field) data visualization »

Membaca citra DICOM (Python)

15Mei10

tampilan file dicom

tampilan file dicom

Melanjutkan cerita tentang membaca file citra berformat DICOM sebelumnya, kali ini mencoba melakukan hal yang sama hanya saja menggunakan python. Sebetulnya sudah ada library untuk tidak hanya membaca melainkan juga mengedit file DICOM yang ditulis murni dengan python (bukan ditulis dalam C++ lalu dibungkus dengan python) yaitu pydicom, tapi bagi saya sekalian melakukan latihan pemrograman python khususnya pembacaan file biner. Selain membaca, saya juga membuat penampil citra dengan menggunakan pygame. Sementara ini masih tampilan per slice sebagai citra grayscale. Pengembangan selanjutnya mau dilanjutkan ke visualisasi 3d dengan menggunakan pyopengl.

Berikut ini kode untuk unit untuk membaca file dicom

#file : dicom.py
import os, sys, struct, string, math

DCMWIDTH = "0028.0010"
DCMHEIGHT = "0028.0011"
DCMALLOCBPP = "0028.0100"
DCMSTOREDBPP = "0028.0101"
DCMSID = "0020.000D"
DCMSTUDYID = "0020.000E"
DCMSPACING = "0028.0030"
DCMPIXELS = "7FE0.0010"
DCMREFUID = "0020.0052"

DCMIMAGEONLY = [DCMWIDTH, DCMHEIGHT, DCMALLOCBPP, DCMSTOREDBPP, DCMSID, DCMSTUDYID, DCMSPACING, DCMPIXELS, DCMREFUID]

class DICOMUnknownChunkError(Exception):
	pass

def readfile(o,cmd):
	return struct.unpack(cmd, o.read(struct.calcsize(cmd)))
	
def load(filename,filter=None):
	fsize = os.path.getsize(filename)
	f = open(filename, "rb")
	f.seek(128)
	cid = f.read(4)
	if cid != "DICM":
		raise Exception("Not a DICOM File")
	data = {}
	while f.tell() < fsize:
		grpelm = readfile(f, 'HH')
		codevr = f.read(2)
		len = readfile(f, 'H')[0]
		id = "%.4x.%.4x" % (grpelm[0], grpelm[1])
		id = id.upper()
		val = None
		if codevr == "UL":
			val = readfile(f, 'I')[0]
		elif codevr == "SL":
			val = readfile(f, 'i')[0]
		elif codevr == "US":
			val = readfile(f, 'H')[0]
		elif codevr == "SS":
			val = readfile(f, 'h')[0]
		elif codevr == "FL":
			val = readfile(f, 'f')[0]
		elif codevr == "FD":
			val = readfile(f, 'd')[0]
		elif codevr in ["LO", "CS", "TM", "DS", "UI", "DA", "PN", "SH", "IS", "ST", "AS", "AE", "AT"]:
			val = str(f.read(len))
		elif codevr in ["OB", "SQ", "UN", "UT"]:
			tmp = readfile(f, 'i')[0]
			val = readfile(f, 'b' * tmp)
		elif codevr == "OW":
			tmp = readfile(f, 'i')[0]
			numelmt = tmp/struct.calcsize('H')
			val = struct.unpack('H'*numelmt, f.read(tmp))
		elif codevr == "OF":
			tmp = readfile(f, 'i')[0]
			numelmt = tmp/struct.calcsize('f')
			val = readfile(f, 'f' * numelmt)
		else:
			raise DICOMUnknownChunkError("unknown chunk %s %s" % (id, codevr))
		if val is not None:
			if filter is None or id in filter:
				data[id] = val
	f.close()
	return data
	
def listseries(dir):
	adir = os.path.abspath(dir)
	result = {}
	for path,dirs,files in os.walk(adir):
		print "looking in", repr(path)
		for file in files:
			try:
				file = os.path.join(path,file)
				data = load(file, [DCMREFUID])
				uid = data[DCMREFUID]
				if uid not in result.keys():
					result[uid] = [file]
				else:
					result[uid].append(file)
			except DICOMUnknownChunkError, e:
				print e
			except:
				pass
	return result

Modul dicom tersebut memiliki dua fungsi utama yaitu load dan listseries. Fungsi load akan menghasilkan objek dictionary yang himpunan kuncinya adalah string dengan format “XXXX.YYYY” yang masing-masing elemennya adalah bilangan heksadesimal. Beberapa telah dibuat konstantanya sehingga tidak perlu mengingat detil pasangan nilainya. Nilai yang ditampung di dictionary ini sesuai dengan nilai yang disimpan di setiap elemen data (angka, string, data biner berbentuk string, dan tuple khusus untuk data piksel yang disimpan sebagai 16-bit integer/word). Fungsi listseries akan menghasilkan objek dictionary yang himpunan kuncinya adalah elemen data Frame of Reference UID, dan tiap elemennya berisi list nama file dari kode tersebut.

modul di atas digunakan pada modul untuk menampilkan menggunakan pygame di kode berikut : 

#file : sliceview.py
from pygame import *
from pygame.locals import *
import numpy as N
import threading
import struct
import dicom, os, sys

def DICOMtosurface(dicomdata):
	w,h = dicomdata[dicom.DCMWIDTH],dicomdata[dicom.DCMHEIGHT]
	
	abpp = dicomdata[dicom.DCMALLOCBPP]
	sbpp = dicomdata[dicom.DCMSTOREDBPP]
	spacing = dicomdata[dicom.DCMSPACING].split("\\")
	pixels = dicomdata[dicom.DCMPIXELS]
	
	#** image.fromstring, fastest (~310ms)
	buffer = ""
	for p in pixels:
		tmp = int(round(p*255.0/65535))
		buffer += struct.pack('BBB',*(tmp,tmp,tmp))
	return image.fromstring(buffer, (w,h), 'RGB')

class DICOMLoader(threading.Thread):
	def __init__(self,filename):
		threading.Thread.__init__(self)
		self.fn = filename
		self.data = None
		self.done = False
	def run(self):
		self.data = dicom.load(self.fn)
		#print self.fn, "done!"
		self.done = True

init()
display.init()
screen = display.set_mode((256, 270), HWSURFACE)
display.update()

imgdicom = []
slices = []
if len(sys.argv)>1:
	slices.append(dicom.load(sys.argv[1], dicom.DCMIMAGEONLY))
	start = time.get_ticks()
	imgdicom.append(DICOMtosurface(slices[0]))
	print time.get_ticks()-start
else:
	#print "use %s <filename>" % sys.argv[0] 
	#exit()
	print "loading default DICOM : "
	threads = []
	for i in range(1,257):
		#hardcode default (hati-hati, edit sebelum digunakan)
		filename = "E:\\data\\dicom\\IM_%.5d.dcm" % i
		loader = DICOMLoader(filename)
		threads.append(loader)
		loader.start()
	alldone = False
	olddone = 0
	while not alldone:
		alldone = True
		numdone = 0
		for t in threads:
			if t.done:
				numdone += 1
			alldone = alldone and t.done
		if numdone > olddone:
			olddone = numdone
			print 100.0*numdone/len(threads),"% loaded"
	for t in threads:
		t.join()
		slices.append(t.data)
		imgdicom.append(None)
	imgdicom[0] = DICOMtosurface(slices[0])

sel = 0
running = True
tfont = font.Font(None, 16)
caption = tfont.render(str(1+sel), 1, (255, 255, 255))

while running:
	screen.fill((127,127,127))
	screen.blit(caption, (0,0))
	screen.blit(imgdicom[sel], (0,caption.get_height()))
	
	display.update()
	for e in event.get():
		if e.type == QUIT or (e.type ==KEYDOWN and e.key == K_ESCAPE):
			running = False
			break
		if e.type == KEYUP:
			if e.key == K_LEFT:
				if sel>0:
					sel -= 1
			elif e.key == K_RIGHT:
				if sel<len(slices)-1:
					sel += 1
			caption = tfont.render(str(1+sel), 1, (255, 255, 255))
			if imgdicom[sel] == None:
				imgdicom[sel] = DICOMtosurface(slices[sel])

Membaca dan mengkonversi data piksel DICOM agar dapat ditampilkan di pygame ternyata tidak secepat dengan delphi. setelah beberapa kali bereksperimen, akhirnya diputuskan bahwa proses pembacaan beberapa file sekaligus dipisah antara proses parsin data DICOM dan konversi ke objek Surface. konversi ke objek surface hanya dilakukan sekali, sedangkan slice lainnya akan dikonversi jika diminta untuk ditampilkan. Proses pembacaan file DICOM pun menggunakan Thread (walaupun sepertinya tidak terlalu signifikan jadinya).

Hal yang sempat menjadi perhatian di kode di atas adalah mentransfer data pixel dari hasil pembacaan citra ke objek Surface dari pygame (pygame.surface.Surface) yang ternyata ada tiga cara. Cara yang paling cepat adalah yang tercantum di kode di atas (tanda #** ). kode di atas juga mengasumsikan data pixel disimpan dalam format 16 bit (unsigned word) dan sudah di-unpack (menggunakan modul struct) menjadi tuple. cara kedua menggunakan modul surfarray yang juga menggunakan modul numpy dan ditunjukkan dalam kode berikut :

def DICOMtosurface(dicomdata):
	w,h = dicomdata[dicom.DCMWIDTH],dicomdata[dicom.DCMHEIGHT]
	
	abpp = dicomdata[dicom.DCMALLOCBPP]
	sbpp = dicomdata[dicom.DCMSTOREDBPP]
	spacing = dicomdata[dicom.DCMSPACING].split("\\")
	pixels = dicomdata[dicom.DCMPIXELS]

	#** use surfarray, slowest (~2100ms)
	surf = surface.Surface((w,h))
	asurf = surfarray.array3d(surf)
	for x, s in enumerate(asurf):
		for y, t in enumerate(s):
			tmp = int(round(pixels[y*w+x]*255.0/65535))
			asurf[x][y] = [tmp, tmp, tmp]
	return surfarray.make_surface(asurf)

Cara terakhir adalah dengan menggunakan metode set_at milik objek Surface seperti diperlihatkan di kode berikut : 

def DICOMtosurface(dicomdata):
	w,h = dicomdata[dicom.DCMWIDTH],dicomdata[dicom.DCMHEIGHT]
	
	abpp = dicomdata[dicom.DCMALLOCBPP]
	sbpp = dicomdata[dicom.DCMSTOREDBPP]
	spacing = dicomdata[dicom.DCMSPACING].split("\\")
	pixels = dicomdata[dicom.DCMPIXELS]

	#** use Surface.set_at slow (390ms)
	surf = surface.Surface((w,h))
	for y in range(h):
		for x in range(w):
			tmp = int(round(pixels[y*w+x]*255.0/65535))
			surf.set_at((x, y), surf.map_rgb((tmp, tmp, tmp, 255)))
	return surf
About these ads

Filed under: programming, riset, work   |  Leave a Comment
Kaitkata:, ,

No Responses Yet to “Membaca citra DICOM (Python)”

Pengumpan untuk Entri ini Trackback Address
  1. Tinggalkan Sebuah Komentar

Tinggalkan Balasan

Fill in your details below or click an icon to log in:

Gravatar
WordPress.com Logo

You are commenting using your WordPress.com account. ( Log Out / Ubah )

Twitter picture

You are commenting using your Twitter account. ( Log Out / Ubah )

Facebook photo

You are commenting using your Facebook account. ( Log Out / Ubah )

Batal

Connecting to %s

Ikuti

Get every new post delivered to your Inbox.

Bergabunglah dengan 1.915 pengikut lainnya.

%d bloggers like this: