Membaca citra DICOM

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!

2 comments

  1. adi · Juni 8, 2010

    salam kenal,

    on this article you wrote :
    ……..namun sepertinya terlalu rumit karena fasilitas yang diberikan terlalu banyak padahal waktu itu kebutuhannya hanya membaca data citranya saja untuk diproses lebih lanjut……..

    kalo boleh tau, waktu itu tujuan pembacaan data citra untuk pengembangan apa ya??

    kalo ngga boleh tau jg gpp😦

    makasih

    • pebbie · Juni 8, 2010

      Salam, waktu itu pake untuk rekonstruksi isosurface jadi saya cuma perlu beberapa field dari DICOMnya saja (spacing, location, seriesID) dan piksel di nilai aslinya (bukan yang sudah jadi bitmap)

Tinggalkan Balasan

Isikan data di bawah atau klik salah satu ikon untuk log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout / Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout / Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout / Ubah )

Foto Google+

You are commenting using your Google+ account. Logout / Ubah )

Connecting to %s