DICOM 이란 Digital Imaging and Communications in Medicine 으로 이 포스팅은 의료영상 이미지 중 CT 이미지 전처리에 대해 다루려고 한다.
.dcm(dicom)파일로 저장된 의료영상 이미지는 float array를 따로 빼내어 딥러닝에 사용한다.
float array를 따로 저장할 수 있는 확장자가 많지만 필자는 주로 tiff 파일로 저장하여 사용한다.
다이콤파일에서 픽셀 어레이만 뽑아서 가져온 후 쌓아서 저장하면 다 아닌가?
라고 생각할 수 있지만, 그게 다가 아니다!
3D영상의 경우, 위처럼 한 폴더 내에 여러개의 .dcm파일이 들어있다. 위 폴더에는 134개의 dcm파일이 들어있음을 볼 수 있다.
그런데, 파일 번호 순서대로 픽셀 어레이를 쌓으면 안된다. 번호 순서대로 쌓으면 위치순서가 엉망진창임을 확인할 수 있다. dcm파일을 읽은 후 InstanceNumber를 뽑아 sort한 후 쌓도록 하자.
from glob import glob
import os
import numpy as np
import pydicom
import skimage.external.tifffile import imsave, imread, imshow
dicom_dir = '/public_data/LIDC-IDRI/LIDC-IDRI-0001/01-01-2000-30178/3000566-03192'
dcm_files = glob(os.path.join(dicom_dir, '*.dcm'))
dcm_files = [pydicom.dcmread(s) for s in dcm_files]
dcm_files.sort(key = lambda x : int(x.InstanceNumber))
dsRef = dcm_files[0]
dcm 파일을 읽으면, 그 안에 많은 key들이 있다. dimension(key : Rows, Columns) 정보 뿐만 아니라, ct vendor사, pixel간 거리(key : PixelSpacing, SliceThickness), pixel array 등이 있으니 필요한 정보가 있으면 key를 뽑아본 후 정보를 확인해 보는 것도 좋은 방법이다.
"""
see all keys in dicom files
"""
print(dsRef)
print(dsRef.ImageType)
for key in dsRef.__dir__():
print(key)
print(key, dsRef.data_element(key))
"""
get dimension and thickness information from the dicom file
"""
dims = (len(dcm_files), int(dsRef.Rows), int(dsRef.Columns))
print('dims(z,x,y) : ', dims)
spacing = (float(dsRef.SliceThickness), float(dsRef.PixelSpacing[0]), float(dsRef.PixelSpacing[1]))
print('thickness(z,x,y) : ', spacing)
이제 dimension정보를 얻었으니 sort한 dicom_files를 토대로 pixel_array값을 불러오도록 하자.
recon_ct = np.zeros(dims, dtype = dsRef.pixel_array.dtype)
for i, df in enumerate(dicom_files):
try :
recon_ct[i,:,:] = df.pixel_array
except :
print(str(i+1).zfill(5) + '.dcm', '**pixel_array_shape Error')
CT영상을 처리하기 위해서는 Hounsfield Units(HU)이 무엇인지 알아야 한다. 이는 X선이 몸을 투과할 때 부위별 흡수정도를 표시한 지표로 CT number라고 부르기도 한다. 물을 0으로 고정하였을 때의 상대적인 흡수량이라고 생각하면 된다(물의 attenuation coefficient에 대한 상대적 비율 * 1000). 우리는 dicom 파일로부터 뽑은 픽셀 어레이를 우선 HU 단위로 정규화 해야한다. 이 때 필요한 것은 dcm파일의 Rescale Slope과 Rescale Intercept 이다.
CT이미지의 경우, 디스크에 저장되는 값과 메모리에 올라오는 값의 표현이 다르게 설정되어 있다. HU는 음수를 포함한 정수값이지만, CT이미지는 일반적으로 unsigned integer인 부호없는 정수로 저장되기 때문이다. 아래의 식은 메모리(output)와 디스크(stored value)에 저장되어 있는 픽셀값의 linear transformation 관계식이다.
(output) = (rescale slope) * (stored value) + (rescale intercept)
우리는 dicom 파일로부터 stored value, rescale slope, rescale intercept을 얻어 위의 식을 계산함으로써 HU에 맞는 값을 얻을 수 있다. (참고 : https://blog.kitware.com/dicom-rescale-intercept-rescale-slope-and-itk/)
"""
Hounsfield Units(HU)
"""
recon_ct = recon_ct.astype(np.int16)
recon_ct[recon_ct == -2000] = 0
intercept = dsRef.RescaleIntercept
slope = dsRef.RescaleSlope
if slope != 1:
recon_ct = slope * recon_ct.astype(np.float64)
recon_ct = recon_ct.astype(np.int16)
recon_ct += np.int16(intercept)
HU에 맞는 값을 얻은 후, 무엇을 보고싶은지에 따라 알맞게 window width와 window level을 조정한다. 예를 들어, -1000~+400까지만 보고싶다면, upper bound = 400, lower bound = -1000으로 설정한 후, 0-1로 normalize하여 사용한다. 이렇게 새로 정의된 recon_ct 를 tiff file로 원하는 path에 맞춰 저장하면 된다.
MIN_BOUND = -1000.0
MAX_BOUND = 400.0
recon_ct = (recon_ct - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
recon_ct[recon_ct>1] = 1.0
recon_ct[recon_ct<0] = 0.0
recon_ct = np.array(recon_ct, dtype = np.float32)
#save array2tiff
target_path = 'D:/data/tiff'
target_name = 'example.tiff'
target_path = os.path.join(target_path, target_name)
imsave(target_path, recon_ct)