wow! post terakhir blog ini waktu musim uts dan sekarang sudah uas. setengah semester gak nulis?
DICOM merupakan standar pengolahan (penyimpanan, komunikasi, penyetakan, dan pemrosesan) untuk keperluan medis. Salah satu subset dari DICOM adalah format berkas untuk citra. Di kegiatan riset KK 2009 kebetulan melibatkan citra DICOM sehingga akhirnya saya ceritakan di sini. Sebetulnya pustaka untuk membaca (dan menulis) citra dalam format DICOM sudah banyak yang tersedia baik komersial maupun yang kodenya tersedia bebas, namun sepertinya terlalu rumit karena fasilitas yang diberikan terlalu banyak padahal waktu itu kebutuhannya hanya membaca data citranya saja untuk diproses lebih lanjut. Akhirnya saya memutuskan untuk mencoba membuat sendiri modul yang khusus untuk membaca citra DICOM.
format penyimpanan citra DICOM sebetulnya cukup sederhana, berkas tersusun atas kumpulan potongan (chunk) informasi yang masing-masing potongan tersebut memiliki metainformasi (header) yang disebut DICOM data element. deskripsinya bisa dilihat pada kode berikut.
type
TDICOMDataElmt = packed record
group, element: word;
VR: array[0..1] of char;
val_length: word;
end;
Deskripsi mengenai maksud potongan informasi tersebut dikodekan dalam bentuk pasangan nilai grup dan elemen (field group dan element) yang bertipe word (16-bit integer) sedangkan field VR merupakan deskripsi tipe data yang dikandung potongan tersebut. field val_length menyatakan panjang sisa informasi potongan tersebut (jika ingin dibaca atau dilewat).
Sebelum membaca urutan potongan informasi seperti yang dijelaskan di atas, ada beberapa hal yang perlu dilakukan yaitu membaca identitas (4 karakter/Bita) pada offset ke-128 (saya juga kurang tahu mengapa 128 Bita pertama dilewat).
Seek( 128, soFromBeginning );
Read( cid[0], 4 );
Kode di bawah merupakan keseluruhan fungsi untuk membaca citra ke dalam tipe internal yaitu TDICOMImage. Bagian yang berisi data piksel biasanya disimpan di VR bertipe “OB”, “OW”, atau “OF”.
//interface
type
TDICOMinfo = record
XSpacing: real;
YSpacing: real;
SeriesID: string;
StudyID: string;
FrameofRefUID: string;
AllocBits, StoredBits: word;
SliceLocation: real;
WindowCenter: real;
WindowWidth: real;
Width, Height: word;
end;
TDICOMimage = record
info: TDICOMInfo;
im: array of array of real;
end;
//...
//implementation
//...
function DICOM_load_image( filename: string; var img: TDICOMimage; skipimage: boolean ): boolean;
var
cid : array[0..3] of char;
stmp : array[0..$FF] of char;
btmp : array of byte;
wtmp : array of word;
ftmp : array of single;
de : TDICOMDataElmt;
svr : string;
itmp : integer;
ustmp : word;
j : integer;
a, b : integer;
isnew : boolean;
pp : TStrings;
begin
result := false;
pp := TStringlist.Create;
try
if not fileexists( filename ) then exit;
isnew := false;
with TFileStream.Create( filename, fmOpenRead ) do begin
Seek( 128, soFromBeginning );
Read( cid[0], 4 );
while Position < Size do begin
Read( de, sizeof( de ) );
svr := string( de.VR );
if ( svr = 'UL' ) or ( svr = 'SL' ) then begin
Read( itmp, sizeof( itmp ) );
end
else
if ( svr = 'US' ) then begin
Read( ustmp, sizeof( ustmp ) );
if de.group = $28 then
if de.element = $10 then begin
if img.info.Width <> ustmp then isnew := true;
img.info.Width := ustmp
end
else
if de.element = $11 then begin
if img.info.Height <> ustmp then isnew := true;
img.info.Height := ustmp
end
else
if de.element = $100 then begin
img.info.AllocBits := ustmp;
end
else
if de.element = $101 then begin
img.info.StoredBits := ustmp;
end;
end
else
if ( svr = 'LO' ) or ( svr = 'CS' ) or ( svr = 'TM' )
or ( svr = 'DS' ) or ( svr = 'UI' ) or ( svr = 'DA' )
or ( svr = 'PN' ) or ( svr = 'SH' ) or ( svr = 'IS' )
or ( svr = 'ST' ) or ( svr = 'AS' ) then begin
Read( stmp[0], de.val_length );
stmp[de.val_length] := #0;
if ( de.group = $20 ) and ( de.element = $D ) then img.info.SeriesID := Strpas( stmp );
if ( de.group = $20 ) and ( de.element = $E ) then img.info.StudyID := Strpas( stmp );
if ( de.group = $20 ) and ( de.element = $52 ) then img.info.FrameofRefUID := Strpas( stmp );
case de.group of
$0020: begin
case de.element of
$1041: begin
img.info.SliceLocation := StrToFloat( strpas( stmp ) );
end;
end;
end;
$0028: begin
case de.element of
$0030: begin
pp.Delimiter := '\';
pp.DelimitedText := strpas( stmp );
img.info.XSpacing := StrToFLoat( pp[0] );
img.info.YSpacing := StrToFloat( pp[1] );
end;
end;
end;
end;
end
else
if ( svr = 'OB' ) or ( svr = 'SQ' ) or ( svr = 'UN' ) or ( svr = 'UT' ) then begin
Read( itmp, sizeof( itmp ) );
Setlength( btmp, itmp );
Read( btmp[0], itmp );
if ( de.group = $7FE0 ) and ( de.element = $10 ) then begin
if not skipimage then begin
if isnew then begin
setlength( img.im, img.info.Height );
for j := 0 to high( img.im ) do
setlength( img.im[j], img.info.Width );
end;
a := 0;
b := 0;
for j := 0 to high( btmp ) do begin
img.im[b, a] := btmp[j];
inc( a );
if a >= img.info.Width then begin
inc( b );
a := 0;
if b >= img.info.Height then break;
end;
end;
end;
end;
end
else
if ( svr = 'OW' ) then begin
Read( itmp, sizeof( itmp ) );
setlength( wtmp, itmp );
Read( wtmp[0], itmp );
if ( de.group = $7FE0 ) and ( de.element = $10 ) then begin
if not skipimage then begin
if isnew then begin
setlength( img.im, img.info.Height );
for j := 0 to high( img.im ) do
setlength( img.im[j], img.info.Width );
end;
a := 0;
b := 0;
for j := 0 to high( wtmp ) do begin
img.im[b, a] := wtmp[j];
inc( a );
if a >= img.info.Width then begin
inc( b );
a := 0;
if b >= img.info.Height then break;
end;
end;
end;
end;
end
else
if ( svr = 'OF' ) then begin
Read( itmp, sizeof( itmp ) );
setlength( ftmp, itmp );
Read( ftmp[0], itmp );
if ( de.group = $7FE0 ) and ( de.element = $10 ) then begin
if not skipimage then begin
if isnew then begin
setlength( img.im, img.info.Height );
for j := 0 to high( img.im ) do
setlength( img.im[j], img.info.Width );
end;
a := 0;
b := 0;
for j := 0 to high( ftmp ) do begin
img.im[b, a] := ftmp[j];
inc( a );
if a >= img.info.Width then begin
inc( b );
a := 0;
if b >= img.info.Height then break;
end;
end;
end;
end;
end
else begin
showmessage( 'unknown ' + svr );
break;
end;
end;
Free;
end;
setlength( btmp, 0 );
setlength( wtmp, 0 );
setlength( ftmp, 0 );
except
exit;
end;
pp.Free;
result := true;
end;
jangan lupa menambahkan unit-unit yang dibutuhkan!
