DICOM의 모든 것 (2)

 

지난 글에는 DICOM의 역사 및 필요성, 확장자에 대한 이해를 주로 다루 었습니다. 

 

이번 글에는 DICOM을 Python으로 불러들이는 방법 및 유용한 라이브러리들을 소개해드리겠습니다.

 

Medical Image 도 다양한 도메인들이 존재합니다. 지난 글에도 말씀드렸다시피 여러 가지의 도메인이나 촬영기법(Sequence)의 이미지를 하나의 확장자로 묶어 사용을 용이하게 하기 위해 DICOM이 나타나게 되었다고 말씀드렸습니다.

 

제가 주로 사용하는 이미지는 자기공명혈관 조영술(Magnetic resonance angiography, MRA)이라는 Sequence 중 3D TOF (time-of-flight) MRA로 촬영된 이미지를 사용하고 있습니다. 해당 촬영기법에 대해서는 추후에 다루어 보도록 하고, 이번 글에서는 Brain MRA 기준으로 Medical Image를 처리하는 방법을 소개해 드리겠습니다.

 

Python으로 DICOM을 Load 하기

 

파이썬으로 DICOM을 load 할 때 사용되는 라이브러리는 Pydicom, SimpleITK 등 다수의 라이브러리가 존재합니다.

 

오늘 사용할 라이브러리는 제가 주로 사용하고 있는 라이브러리인 Pydicom을 소개해 드리려고 합니다.

 

Pydicom을 설치하는 방법은 따로 해당 라이브러리 Docs에서 확인하시는 것이 편할듯 합니다. 

 

 

먼저 한장의 이미지를 Load 하는 법을 살표 보겠습니다.

 

Pydicom으로 Dicom 파일 Load 하기

import pydicom

path = "본인의 DICOM 데이터 path/***.dcm"
slice = pydicom.dcmread(path)

print(slice)

 

해당 코드를 컴파일 하게되면 아래와 같은 결과 값을 얻을 수 있습니다. (내용은 개인 정보를 익명화 하였지만 혹시 모를 상황을 위해 블러 처리하였습니다.)

 

Code 캄파일 결과

Dataset file_meta 정보들을 얻게 되며 다양한 Tag 값들을 확인할 수 있습니다.

 

그러면 이제 불러온 데이터 중 이미지 픽셀 데이터에 접근하여 시각화해 보겠습니다.

 

Load 된 데이터 시각화 하기

import pydicom
import matplotlib.pyplot as plt

path = "본인의 DICOM 데이터 path/***.dcm"
data_slice = pydicom.dcmread(path)

plt.imshow(data_slice.pixel_array, cmap='gray')
plt.show()

 

불러온 데이터를 저장한 객체인 data_slice 변수에 pixel_array에 접근하여 시각화를 진행해 보았습니다.

 

Matplot의 Imshow함수의 Color Map Default 값은 "viridis"이기 때문에 cmap를 설정하지 않고 시각화를 진행하면 왼쪽과 같은 이미지가 출력되는 것을 확인할 수 있습니다.

 

따라서 일반적으로 DICOM Viewer에서 보는 색감으로 보고 싶은 분께서는 cmap =  'gray'을 사용하시면 될 것 같습니다.

 

 

 

이번에는 DICOM을 여러 장 불러오는 방법에 대해서도 진행해 보겠습니다.

 

Brain MRA와 같은 경우는 한 장 한 장의 사진으로 병변을 확인하기는 어려움이 있습니다.

 

따라서 이를 Post Processing를 통해 MIP(Maximum intensity projection)를 통해 뇌졸중(뇌 동맥류,,, etc)의 위치 등을 확인합니다. 그러기 위해선 여러 장으로 구성된 뇌 구조를 하나의 리스트에 모을 필요성이 있습니다.

Pydicom으로 다수의 DICOM파일 Load

import pydicom
import matplotlib.pyplot as plt
import os
import glob
from tqdm import tqdm


path = "다수의 DICOM이 존재하는 폴더 path"
datasets = [pydicom.dcmread(slice) for slice in tqdm(glob.glob(os.path.join(path, '*.dcm')))]

 

하나의 리스트 변수에 다수의 DICOM 파일들을 Load 해보았습니다. 이때 사용한 라이브러리들 중 os, glob 폴더의 파일들 중 확장자가 DICOM이 아닌 파일들이 있을 수 있기 때문에 "glob.glob(os.path.join(path, '*. dcm'))"으로 확장자가 ". dcm"만 불러올 수 있도록 하였고, tqdm은 파일을 불러올 때 진행률을 보기 위한 Progress bar을 설정하기 위한 라이브러리입니다. 자세한 라이브러리 내용은 따로 찾아보시면 유용할 것 같다고 생각합니다. 

 

이제 하나의 리스트로 다수의 DICOM 파일들을 불러왔으니 이걸 시각화를 한번 진행해보겠습니다.

Anatomical Plane

진행하기에 앞서 해부학적 평면에 대한 내용을 조금 다루고 지나가겠습니다.

 

아래의 그림은 해부학적 위치에 대한 그림을 나타냅니다.

 

 

Anatomical Position

 

An anatomical plane is a hypothetical plane used to transect the body, in order to describe the location of structures or the direction of movements. In human and animal anatomy, three principal planes are used:

• The sagittal plane or lateral plane (longitudinal, anteroposterior) is a plane parallel to the sagittal suture. It divides the body into left and right.
• The coronal plane or frontal plane (vertical) divides the body into dorsal and ventral (back and front, or posterior and anterior) portions.
• The transverse plane or axial plane (horizontal) divides the body into cranial and caudal (head and tail) portions.

 

해부학적 평면이라는 내용의 위키 내용을 첨부하였습니다. 번역보다는 원어 그대로 보는 것이 이해가 편할 것 같아 원어 그대로 올립니다. (사실 영어 잘못해서 ㅎㅎ)

 

 

Multi Planar Reconstruction(MPR)

해부학적 평면을 기반으로 하여 MR, CT 등의 시퀀스 이미지를 시각화하는 방법인 Multi Planar Reconstruction을 통해 시각화를 진행해 보겠습니다.

 

아래의 그림은 Radiant라는 프로그램을 통해 MPR기반의 시각화 기능을 나타냅니다.

위와 같이 Python에서 MPR 기반 시각화를 진행해 보겠습니다. 

 

Pydicom으로 MPR 기반 시각화 하기

 

import imp
import pydicom
import matplotlib.pyplot as plt
import os
import glob
from tqdm import tqdm
import numpy as np


path = "다수의 DICOM이 존재하는 폴더 path"
files = [pydicom.dcmread(slice) for slice in tqdm(glob.glob(os.path.join(path, '*.dcm')))]


# skip files with no SliceLocation (eg scout views)
slices = []
skipcount = 0
for f in files:
    if hasattr(f, 'SliceLocation'):
        slices.append(f)
    else:
        skipcount = skipcount + 1

print("skipped, no SliceLocation: {}".format(skipcount))

# ensure they are in the correct order
slices = sorted(slices, key=lambda s: s.SliceLocation)

# pixel aspects, assuming all slices are the same
ps = slices[0].PixelSpacing
ss = slices[0].SliceThickness
ax_aspect = ps[1]/ps[0]
sag_aspect = ps[1]/ss
cor_aspect = ss/ps[0]

# create 3D array
img_shape = list(slices[0].pixel_array.shape)
img_shape.append(len(slices))
img3d = np.zeros(img_shape)

# fill 3D array with the images from the files
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:, :, i] = img2d

# plot 3 orthogonal slices
a1 = plt.subplot(2, 2, 1)
plt.imshow(img3d[:, :, img_shape[2]//2],cmap='gray')
a1.set_aspect(ax_aspect)

a2 = plt.subplot(2, 2, 2)
plt.imshow(img3d[:, img_shape[1]//2, :],cmap='gray')
a2.set_aspect(sag_aspect)

a3 = plt.subplot(2, 2, 3)
plt.imshow(img3d[img_shape[0]//2, :, :].T,cmap='gray')
a3.set_aspect(cor_aspect)

plt.show()

 

해당 코드는 pydicom Docs에서 퍼온 내용입니다. 다수의 DICOM을 Load 한 변수를 사용하여 "Slice Location" Tag 유무를 확인하고 해당 태그의 Value를 통해 DICOM을 정렬하는 모습을 볼 수 있습니다.

# skip files with no SliceLocation (eg scout views)
slices = []
skipcount = 0
for f in files:
    if hasattr(f, 'SliceLocation'):
        slices.append(f)
    else:
        skipcount = skipcount + 1

print("skipped, no SliceLocation: {}".format(skipcount))

 

해당 코드는 Axial, Sagittal, Cornal 축을 설정하는 코드입니다.

# pixel aspects, assuming all slices are the same
ps = slices[0].PixelSpacing
ss = slices[0].SliceThickness
ax_aspect = ps[1]/ps[0]
sag_aspect = ps[1]/ss
cor_aspect = ss/ps[0]

다음으로 위치한 코드는 3차원 Array를 생성하고 기존의 2D slice로 존재하는 데이터를 3D Voxel로 재구성하는 모습을 볼 수 있습니다. Voxel의 개념은 추후에 따로 다루어보겠습니다.

# create 3D array
img_shape = list(slices[0].pixel_array.shape)
img_shape.append(len(slices))
img3d = np.zeros(img_shape)

# fill 3D array with the images from the files
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:, :, i] = img2d

마지막은 3차원으로 구성한 데이터를 시각화하는 부분입니다. Subplot를 통해 하나의 창에 다수의 이미지를 시각화할 수 있게 하였고 Img_shape [0], Img_shape [1], Img_shape [2]를 통해 각각의 해부학적 평면의 중간에 위치한 이미지를 시각화할 수 있도록 하였습니다.

# plot 3 orthogonal slices
a1 = plt.subplot(2, 2, 1)
plt.imshow(img3d[:, :, img_shape[2]//2],cmap='gray')
a1.set_aspect(ax_aspect)

a2 = plt.subplot(2, 2, 2)
plt.imshow(img3d[:, img_shape[1]//2, :],cmap='gray')
a2.set_aspect(sag_aspect)

a3 = plt.subplot(2, 2, 3)
plt.imshow(img3d[img_shape[0]//2, :, :].T,cmap='gray')
a3.set_aspect(cor_aspect)

plt.show()

다음은 전체 코드의 결과 그림을 나타냅니다.

pydicom을 통한 MPR기반 시각화 결과

결과물을 자세히 보면 Radiant보다 늘어나 있거나 회전되어 있는 모습을 볼 수 있습니다. 해당 문제는 다음 글에서 다루어 보겠습니다.

 

이번 글은 Pydicom을 통해 DICOM을 다루는 방법에 대해 설명하여 보았습니다.

 

긴 글 읽어 주셔서 감사합니다.