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
GAIBlog 