[SimpleITK] Python으로 medical image processing하기 - (1) Load, Resample

의료영상의 전처리에 정말 많이 쓰이는 library인 것 같은데, 정보가 많이 없어 힘들었다.

SimpleITK library를 이용해 Python으로 이미지를 resample하는 법과, 이후에 가능하다면 registration 하는 법까지 정리해 보려고 한다.

 

SimpleITK를 이용해 Python으로 processing 하기

[SimpleITK] Python으로 medical image processing하기 - (1) Load, Resample

[SimpleITK] Python으로 medical image processing하기 - (2) Registration 1

[SimpleITK] Python으로 medical image processing하기 - (3) Registration 2

[SimpleITK] Python으로 medical image processing하기 - (4) Visualization

 

 

install & import SimpleITK

pip install SimpleITK
pip install SimpleITK-SimpleElastix

import SimpleITK as sitk

 

 

numpy array를 sitk.Image 객체로 바꾸기

sitk.GetImageFromArray() method를 이용해 numpy array를 sitk.Image 객체로 바꿔줄 수 있다.

 

- Set~ method를 이용해 다양한 attribute(origin, spacing 등)을 설정해 줄 수 있으며,

- Get~ method를 이용해 객체의 attribute를 확인할 수 있다.

 

- 전체 method와 예시(링크)

    - ex) sitk_volume.GetOrigin() : sitk_volume의 origin을 return한다.

             sitk_volume.GetSize() : sitk_volume의 size를 return한다.

             sitk_volume.GetDepth() : sitk_volume의 depth(마지막 dim)를 return한다.

             sitk_volume.GetPixel(pixel_location) : pixel_location 위치의 pixel 값을 return한다.

             sitk_volume.SetOrigin(origin) : sitk_volume의 origin을 설정한다.

def convert_to_sitk(img_volume, spacing, origin, direction=None):
    sitk_volume = sitk.GetImageFromArray(img_volume)
    sitk_volume.SetOrigin(origin)
    sitk_volume.SetSpacing(spacing)
    if direction:
    	sitk_volume.SetDirection(direction)
    return sitk_volume

 

GetImageFromArray method는 input으로 들어오는 numpy image를 [z, y, x] 순서로 생각한다.

즉, input image가 (2, 1, 0) 순서로 transpose 된다.

import numpy as np
import SimpleITK as sitk

x = np.random.rand(3, 4, 5)
x.shape #(3, 4, 5)

x_sitk = sitk.GetImageFromArray(x)
x_sitk.GetSize() #(5, 4, 3)

 

이를 원치 않으면 GetImageFromArray method를 호출하기 전에 img_volume을 (2, 1, 0) 순서로 transpose 해주면 된다.

x = np.random.rand(3, 4, 5)
x.shape #(3, 4, 5)

x = x.transpose(2, 1, 0)
x_sitk = sitk.GetImageFromArray(x)
x_sitk.GetSize() #(3, 4, 5)

 

 

 

sitk.Image 객체를 numpy array로 바꾸기

sitk.Image 객체를 다시 numpy array로 바꿀 때는 sitk.GetArrayFromImage() method를 사용하면 된다.

img_volume = sitk.GetArrayFromImage(sitk_volume)

 

마찬가지로 dimension이 (2, 1, 0) 순서로 transpose 되어 나온다.

numpy array에서 sitk.Image 객체로 바꿔주기 전 transpose를 했다면,

sitk.Image 객체에서 numpy array로 바꿔줄 때 역시 transpose를 해 주어야 한다.

 

나는 다음과 같이 둘을 전환하는 함수를 만들어 사용했다.

def convert_to_sitk(img_volume, spacing, origin, direction=None):
    """convert numpy volume to sitk image"""
    img_volume = img_volume.transpose(2, 1, 0) #since sitk assume numpy in [z, y, x]
    sitk_volume = sitk.GetImageFromArray(volume.astype(np.float64))
    sitk_volume.SetOrigin(origin)
    sitk_volume.SetSpacing(spacing)
    if direction:
        sitk_volume.SetDirection(direction)
    return sitk_volume

def convert_to_numpy(sitk_volume):
    """convert sitk image to numpy volume"""
    img_volume = sitk.GetArrayFromImage(sitk_volume)
    img_volume = img_volume.transpose(2, 1, 0)
    return img_volume

 

 

 

nifti/dicom 파일을 바로 불러오기

numpy array가 아닌 .nii.gz / .dcm 파일을 바로 sitk.Image 객체로 바꿀 수도 있다.

 

nifti 파일 불러오기

nii_sitk = sitk.ReadImage('image.nii.gz')

 

dicom series 불러오기

dicom_path = './dicom_path'

series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(dicom_path)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(dicom_path, series_ids[0])

series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)
dicom_sitk = series_reader.Execute()

 

단일 dicom 파일 불러오기

dcm_image = sitk.ReadImage('image.dcm')

 

 

Resample

- SimpleITK는 다음과 같은 순서로 전처리를 수행할 수 있다.

 

 1. 전처리를 수행할 객체를 생성

 

 2.  객체에 parameter 설정

        - 이때 일부 parameter들은(direction 등) 설정하지 않으면 default 값을 사용한다.

        - resample에 사용할 parameter들은 대부분 dicom header에서 가져올 수 있다. 이 글의 마지막 파트를 참고

        - default_value는 background pixel intensity value이다. (ex: CT image의 경우 -1024, MR image의 경우 0으로 설정)

 

3. Execute() method를 이용해 전처리 수행

def resample(sitk_volume, new_spacing, new_size, default_value=0):
    """1) Create resampler"""
    resample = sitk.ResampleImageFilter() 
    
    """2) Set parameters"""
    #set interpolation method, output direction, default pixel value
    resample.SetInterpolator(sitk.sitkLinear)
    resample.SetOutputDirection(sitk_volume.GetDirection())
    resample.SetDefaultPixelValue(default_value)
    
    #set output spacing
    new_spacing = np.array(new_spacing)
    resample.SetOutputSpacing(new_spacing)
    
    #set output size and origin
    old_size = np.array(sitk_volume.GetSize())
    old_spacing = np.array(sitk_volume.GetSpacing())
    new_size_no_shift = np.int16(np.ceil(old_size*old_spacing/new_spacing))
    old_origin = np.array(sitk_volume.GetOrigin())
    
    shift_amount = np.int16(np.floor((new_size_no_shift - new_size)/2))*new_spacing
    new_origin = old_origin + shift_amount
    
    new_size = [int(s) for s in new_size]
    resample.SetSize(new_size)
    resample.SetOutputOrigin(new_origin)
    
    """3) execute"""
    new_volume = resample.Execute(sitk_volume)
    return new_volume

 

spacing, origin 등에 대한 설명은 다음 그림을 참고

https://simpleelastix.readthedocs.io/Introduction.html

 

 

---

 

참고 : official docs (https://simpleitk.org/SPIE2018_COURSE/images_and_resampling.pdf)