醫學圖像處理流程

Tips: 最近在做醫學圖像預處理（CT/PET），涉及到了一些盲點和知識點，在這做一些總結。

一、數據格式

1. DICOM

   DICOM是醫學圖像中的標准文件，包含了許多元數據信息，這些信息具體可以分為以下四類：

   • Patient

   • Study

   • Series

   • Image

   每一個DICOM Tag都是由兩個十六進制數的組合來確定的，分別為Group和Element。如(0010,0010)這個Tag表示的是Patient’s Name，它存儲着這張DICOM圖像的患者姓名。

   

   每個病人的每個模態的有幾十到幾百的dcm數據文件（slices）

2. mhd

   mhd格式是另外一種數據格式。每個病人一個mhd文件和一個同名的raw文件，一個mhd通常有幾百兆，對應的raw文件只有1kb。mhd文件需要借助python的SimpleITK包來處理。SimpleITK示例代碼如下:

   1. import SimpleITK as sitk
   2. itk_img = sitk.ReadImage(img_file)
   3. img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
   4. num_z, height, width = img_array.shape #heightXwidth constitute the transverse plane
   5. origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm)
   6. spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)
   

   需要注意的是，SimpleITK的img_array的數組不是直接的像素值，而是相對於CT掃描中原點位置的差值，需要做進一步轉換

3. NIFIT

   醫學影像早期使用的是DICOM標准，基本上各家廠商都會使用符合DICOM標准的產品，但是這個標准對於數據分析並不方便。在神經影像興起時就誕生了各種各樣的數據存儲標准，比如analyze。NIFTI出現的原因是原來一種圖像格式是ANALYZE 7.5 format，但是這個圖像格式缺少一些信息，比如沒有方向信息，病人的左右方位等，如果需要包括額外的信息，就需要一個額外的文件，比如ANALYZE7.5就需要一對<.hdr, .img>文件來保存圖像的完整信息。因此，解決這個問題Data Format Working Group (DFWG) 將圖像格式完整的定義為NIFTI（Neuroimaging Informatics Technology Initiative）格式。放射學和神經學醫生使用軟件和使用習慣會有一些不同，所以為了統一格式，DFWG提出NIFTI格式圖像。

   NIFTI格式的讀取軟件，也應該可以讀取<.hdr, .img>文件對。

   NIFTI格式中保存的信息是，一共可以保存7維數據。第1,2,3維度（0維統計使用維度個數）保留給空間維度x，y，z，而第四維度留給時間維度t。第5,6,7維度可供用戶其他使用，但是第五維度也有一些預定義的用途，例如存儲特定體素分布的參數或存儲矢量數據。

二、圖像處理

1. dicom數據讀取

   def readdcm(filepath):
      #filepath = "./T2"
       series_id = sitk.ImageSeriesReader.GetGDCMSeriesIDs(filepath)
       series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(filepath, series_id[0])
       series_reader = sitk.ImageSeriesReader() #讀# 取數據端口
       series_reader.SetFileNames(series_file_names)  #讀取名稱
       images = series_reader.Execute()#讀取數據
       #print(images.GetSpacing())
       #sitk.WriteImage(images, "T2_1.nii.gz")#保存為nii
       return images
   
   

2. 其他格式

   image = sitk.ReadImage(filepath)
   

三、圖像預處理

對於不同的數據類型重采樣的方法和目的都不相同。例如在遙感中，重采樣是從高分辨率遙感影像中提取出低分辨率影像的過程；在數據挖掘中，重采樣是指為了解決訓練數據類別不均衡，通過在訓練期間通過增加小樣本的數量或者減少大樣本的數量保持樣本類別均衡的算法；在醫療圖像中，重采樣是指將醫療圖像中大小不同的體素歸一化到相同的大小。體素是體積元素（Volume Pixel）的簡稱，一張3D醫療圖像可以看成是由若干個體素構成的

不同掃描面的像素尺寸、粗細粒度是不同的。這不利於我們進行CNN任務，我們可以使用同構采樣。

一個掃描面的像素區間可能是[2.5,0.5,0.5],即切片之間的距離為2.5mm。可能另外一個掃描面的范圍是[1.5,0.725,0.725]。這可能不利於自動分析。

常見的處理方法是從全數據集中以固定的同構分辨率重新采樣，將所有的東西采樣為1mmx1mmx1mm像素。

對於醫療圖像重采樣，可以分成三個步驟：首先使用SimpleITK讀取數據，獲得原始圖像的對應的Spacing以及Size，得到圖像原始的大小；然后圖像原始的大小除以新Spacing得到新Size；最后將新Spacing得到新Size賦值到讀取的數據即可。

對於一張大小為128128的彩色圖像，其在計算機中可以表示為1281283的矩陣，其中每一個像素點的取值范圍為0-255，不同的數值代表不同的亮度。但是對於醫療圖像其是由若干個slice組成的，假設每一個slice的大小為512512的單通道的圖像，其中每一個像素點表示的是一個體素的取值，其范圍可以-1000~2000之間。接下來通過以胰腺分割數據集中PANCREAS_0015.nii.gz為例，對醫療圖像中體素這個概念進行講解。Spacing(0.78125, 0.78125, 1.0)表示的是原始圖像體素的大小，也可以將Spacing想象成大小為(0.78125, 0.78125, 1.0)的長方體。而原始圖像的Size為 (512, 512, 247)，表示的是原始在X軸，Y軸，Z軸中體素的個數。原始圖像的size對應的Spacing既可以得到真實3D圖像大小（5120.78125，5120.78125，2471 ），在圖像重采樣只是修改體素的大小，而真實3D圖像大小是保持不變的，因此假設我們將Spacing修改成(1.0, 1.0, 2.0)的時候，則修改之后其對應的size應該為（（5120.78125）/ 1.0，（5120.78125）/ 1.0，（247*1 ））即（400, 400, 124）。

def transform(image,newSpacing, resamplemethod=sitk.sitkNearestNeighbor):
    # 設置一個Filter
    resample = sitk.ResampleImageFilter()
    # 原來的體素塊尺寸
    originSize = image.GetSize()
    # 原來的體素間的距離
    originSpacing = image.GetSpacing()
    #print(originSpacing)
    # newSize = np.array(newSize, float)
    # newSpacing = np.array(newSpacing, float)
    # factor = originSpacing / newSpacing
    # newSize = originSize * factor
    # newSpacing = newSpacing.astype(np.int)
    newSize = [
        int(np.round(originSize[0] * originSpacing[0] / newSpacing[0])),
        int(np.round(originSize[1] * originSpacing[1] / newSpacing[1])),
        int(np.round(originSize[2] * originSpacing[2] / newSpacing[2]))
    ]

    # 默認像素值（2）
    # resample.SetDefaultPixelValue(image.GetPi);
    # 沿着x,y,z,的spacing（3）
    # The sampling grid of the output space is specified with the spacing along each dimension and the origin.
    resample.SetOutputSpacing(newSpacing)
    # 設置original
    resample.SetOutputOrigin(image.GetOrigin())
    # 設置方向
    resample.SetOutputDirection(image.GetDirection())
    resample.SetSize(newSize)
    # 設置插值方式
    resample.SetInterpolator(resamplemethod)
    # 設置transform
    resample.SetTransform(sitk.Euler3DTransform())
    # 默認像素值   resample.SetDefaultPixelValue(image.GetPixelIDValue())
    return resample.Execute(image)

if __name__ == '__main__':
    for i in range(1, 103):
        # Dicom序列所在文件夾路徑
        strr = str(i).zfill(3)
        print("Start_", i, '_converse')
        filepath = ''
        filepathLabel = ''

        image = sitk.ReadImage(filepath)  # 讀取image文件
        # print(image.GetOrigin())
        # print(image.GetSpacing())
        # print(image.GetSize())
        label=sitk.ReadImage(filepathLabel) #讀取mask

        newImage=transform(image,[1,1,1],sitk.sitkLinear)
        newLabel=transform(label,[1,1,1])# 不添加插值方式時，默認為最近鄰插值
        # mask請用最近鄰插值，image用線性插值
        name_data = ''# image存儲path
        name_label='' # label存儲path

        sitk.WriteImage(newImage, name_data)
        sitk.WriteImage(newLabel, name_label)