隨著最近幾大肺部圖像處理相關(guān)的競賽的推出,，如LUNA16,、Kaggle Data Science Bowl，AI領(lǐng)域的科研人員對(duì)肺部CT圖像變得越來越熟悉,，尤其是DICOM序列,，以及這些競賽官方所提供的mhd數(shù)據(jù)格式。
ITK是一個(gè)功能很強(qiáng)大的醫(yī)學(xué)圖像處理公開庫,，搭配VTK用以顯示圖像,，可以實(shí)現(xiàn)幾乎所有醫(yī)學(xué)圖像處理的功能需要。ITK通常以C++包進(jìn)行提供,，當(dāng)然也可以自己編譯為Python包,，不過編譯過程比較繁瑣耗時(shí)，而且很容易踩坑,。但I(xiàn)TK官方進(jìn)行的Python封裝SimpleITK,，則直接可以拿來使用；雖然有部分ITK的功能沒有包含,，但已基本夠用了,。我們?cè)谔幚磲t(yī)學(xué)圖像時(shí)，使用的基本都是SimpleITK,。
本文就簡單總結(jié)一下我們?cè)谔幚磉@些圖像時(shí)的經(jīng)驗(yàn),，以便備忘，并為后來者參考,。
1.讀取文件
讀取DICOM序列
醫(yī)學(xué)圖像中一個(gè)CT序列包含很多張圖片,，即一個(gè)case包含許多slice,，使用SimpleITK可以直接讀取一個(gè)序列，并方便地得到各種參數(shù),，將圖像數(shù)據(jù)轉(zhuǎn)換成numpy Array：
import SimpleITK as sitk
import numpy as np
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames(case_path)
reader.SetFileNames(dicom_names)
image = reader.Execute()
image_array = sitk.GetArrayFromImage(image) # z, y, x
origin = image.GetOrigin() # x, y, z
spacing = image.GetSpacing() # x, y, z
需要注意的是,，SimpleITK讀取的圖像數(shù)據(jù)的坐標(biāo)順序?yàn)閦yx，即從多少張切片到單張切片的寬和高,；而據(jù)SimpleITK Image獲取的origin和spacing的坐標(biāo)順序則是xyz,。這些需要特別注意。
讀取DICOM單張圖片
可以將一個(gè)DICOM序列作為一個(gè)整體一次讀入,，也可以一張一張地讀入每張切片：
import SimpleITK as sitk
import numpy as np
image = sitk.ReadImage(slice_path)
image_array = sitk.GetArrayFromImage(image) # z, y, x
這里需要注意的除了坐標(biāo)順序是zyx之外,，還需注意，即使讀取單張切片,，所得到的結(jié)果也是3維的,，只不過第一維是1。
讀取mhd文件
涉及DICOM序列時(shí),，為了傳輸方便（從上百個(gè)dcm文件到一個(gè)mhd文件）,，很多情況下以mhd文件格式進(jìn)行呈現(xiàn)，不過mhd文件只是一個(gè)很小的包含數(shù)據(jù)信息的文件,，同時(shí)搭配的通常還有一個(gè)二進(jìn)制的數(shù)據(jù)文件（格式為raw或zraw等）,。使用SimpleITK讀取這種文件也比較方便。
import SimpleITK as sitk
import numpy as np
image = sitk.ReadImage(mhd_path)
image_array = sitk.GetArrayFromImage(image) # z, y, x
origin = image.GetOrigin() # x, y, z
...
有時(shí)候不想整個(gè)讀取數(shù)據(jù)（因?yàn)楸容^大,，讀取處理比較慢）,，想要讀取的只是一些基本信息，比如origin,、spacing等等,。這時(shí)可以只讀取mhd文件，據(jù)此獲取信息,，讀取方法比較簡單,，不再贅述。
2.處理文件
處理DICOM文件主要有插值等操作,，可以直接使用SimpleITK（或者說是ITK）的相關(guān)函數(shù),，并通過pipeline結(jié)構(gòu)進(jìn)行處理。
插值
import SimpleITK as sitk
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames(case_path) reader.SetFileNames(dicom_names)
image = reader.Execute()
resample = sitk.ResampleImageFilter()
resample.SetOutputDirection(image.GetDirection())
resample.SetOutputOrigin(image.GetOrigin())
newspacing = [1, 1, 1]
resample.SetOutputSpacing(newspacing)
newimage = resample.Execute(image)
