使用Python對Dicom文件進行讀取與寫入

Pydicom

單張影像的讀取

使用 pydicom.dcmread() 函數進行單張影像的讀取,返回一個pydicom.dataset.FileDataset對象.html

import os
import pydicom
# 調用本地的 dicom file 
folder_path = r"D:\Files\Data\Materials"
file_name = "PA1_0001.dcm"
file_path = os.path.join(folder_path,file_name)
ds = pydicom.dcmread(file_path)

在一些特殊狀況下,好比直接讀取從醫院拿到的數據(未經任何處理)時,可能會發生如下報錯:python

raise InvalidDicomError("File is missing DICOM File Meta Information "
pydicom.errors.InvalidDicomError: File is missing DICOM File Meta Information header or the 'DICM' prefix is missing from the header. Use force=True to force reading.

能夠看到,因爲缺失文件元信息頭,沒法直接讀取,只能強行讀取.這種狀況能夠直接根據提示,調整命令爲:git

ds = pydicom.dcmread(file_path,force=True)

但後續還會碰到:github

AttributeError: 'Dataset' object has no attribute 'TransferSyntaxUID'

在網上檢索後發現,能夠經過設置TransferSyntaxUID來解決問題:markdown

ds.file_meta.TransferSyntaxUID = pydicom.uid.ImplicitVRLittleEndian

這樣就大功告成了(這裏實際上就提早接觸到了下面讀取Dicom Tags的內容了)函數

一些簡單處理

讀取成功後,咱們能夠對 Dicom文件 進行一些簡單的處理工具

讀取並編輯Dicom Tags

能夠經過兩種方法來讀取Tag的值ui

  1. 使用的Tag的Description
print(ds.PatientID,ds.StudyDate,ds.Modality)
  1. 使用 ds.get() 函數. 函數內參數採用的是Tag ID.幾種簡單的打開Dicom文件的軟件(如RadiAnt DICOM Viewer)均可以直接看到.這裏再也不贅述.
ds.get(0x00100020) # 這裏獲得的是PatientID 

讀取到相應的Tag值後, 也能夠將其餘的值寫入這些Tag.只要最後保存一下就能夠了.spa

藉助Numpy與PIL.Image

讀取Dicom文件後,能夠藉助Numpy以及圖像處理庫(如PIL.Image)來進行簡單的處理..net

  1. 藉助Numpy
import numpy as np
data = np.array(ds.pixel_array)

注意這裏使用的是 np.array() 而不是 np.asarray(). 由於前者的更改並不會帶來原pixel_array的改變.
在轉化爲ndarray後 能夠直接進行簡單的切割和鏈接,好比截取某一部分和將兩張圖像拼在一塊兒等,以後再寫入並保存下來便可.

  1. 藉助PIL.Image
from PIL import Image
data_img = Image.fromarray(ds.pixel_array)
data_img_rotated = data_img.rotate(angle=45,resample=Image.BICUBIC,fillcolor=data_img.getpixel((0,0)))

這裏展現的是旋轉, 還有其餘功能如resize等.
須要注意的是,從Numpy的ndarray轉化爲Image時,通常會發生變化:

print(data.dtype)  # int16
data_rotated = np.array(data_img_rotated)
print(data_img)   # int32

只須要指定參數就能夠解決了

data_rotated = np.array(data_img_rotated,dtype = np.int16)

可視化

簡單的可視化Pydicom沒有直接的實現方法,咱們能夠經過上面藉助Matplotlib以及Image模塊來實現.但效果有限.

  1. 藉助 Matplotlib (Pydicom官方文檔中使用的方法)
from matplotlib import pyplot
pyplot.imshow(ds.pixel_array,cmap=pyplot.cm.bone)
pyplot.show()

效果如圖所示:
使用Matplotlib進行Dicom文件的可視化
但真實的圖像是:

真實MR圖像
顯然顏色是有區別的.致使這種差異的緣由是pyplot函數使用的cm也就是"color map" 是簡單的"bone" 並不能知足醫學圖像的要求.

  1. 藉助Image模塊
data_img.show()

一條指令便可,可是效果不好,如圖所示:

使用Image模塊進行可視化
綜合來看,兩種方法都不是很好.

單張影像的寫入

通過上面對Tag值的修改, 對圖像的切割, 旋轉等操做.最後須要從新寫入該Dicom文件.

ds.PixelData = data_rotated.tobytes()
ds.Rows,ds.Columns = data_rotated.shape
new_name = "dicom_rotated.dcm"
ds.save_as(os.path.join(folder_path,new_name))

SimpleITK

SimpleITK 是從基於C++的ITK遷移到Python的,因此不少方法的使用都跟C++很類似.

import SimpleITK as sitk

單張影像的讀取

有兩種方法:

  1. sitk.ReadImage()
    這種方法直接返回image對象,簡單易懂.可是沒法讀取Tag的值.
img = sitk.ReadImage(file_path)
print(type(img)) # <class 'SimpleITK.SimpleITK.Image'>
  1. sitk.ImageFileReader()
    這種方法比較像C++的操做風格,須要先初始化一個對象,而後設置一些參數,最後返回image.相對更復雜,但能夠操做的點比較多
file_reader = sitk.ImageFileReader()
file_reader.SetFileName(file_path) #這裏只顯示了必需的,還有不少能夠設置的參數
data = file_reader.Execute()
# 使用這種方法讀取Dicom的Tag Value
for key in file_reader.GetMetaDataKeys():
    print(key,file_reader.GetMetaData(key))

以上兩種方法返回的都是三維的對象,這與Pydicom有很大的不一樣.

data_np = sitk.GetArrayFromImage(data)
print(data_np.shape)  # (1, 512, 512) = (Slice index, Rows, Columns)

序列讀取

序列讀取的方法與單張圖像讀取的第二種方法很類似.
(暫且只發現了一種方法讀取序列,若是還有其餘方法,請在評論區予以補充,感謝!)

series_reader = sitk.ImageSeriesReader()
fileNames = series_reader.GetGDCMSeriesFileNames(folder_name)
series_reader.SetFileNames(fileNames)
images = series_reader.Execute()

一樣,返回的也是三維的對象.

一些簡單操做

SimpleITK 包含不少圖像處理如濾波的工具,這裏簡單介紹一個邊緣檢測工具和可視化工具

邊緣檢測

以Canny邊緣檢測算子爲例,與讀取單張圖像相似,一樣有兩種方式:

  1. sitk.CannyEdgeDetection()
    因爲濾波的對象必須是32位圖像或者其餘格式, 須要經過 sitk.Cast() 轉換. 以後能夠再轉換回原格式.
data_32 = sitk.Cast(data,sitk.sitkFloat32)
data_edge_1 = sitk.CannyEdgeDetection(data_32,5,30,[5]*3,[0.8]*3)
  1. sitk.CannyEdgeDetectionImageFilter()
    這個操做相對麻煩一些
Canny = sitk.CannyEdgeDetectionImageFilter()
Canny.SetLowerThreshold(5)
Canny.SetUpperThreshold(30)
Canny.SetVariance([5]*3)
Canny.SetMaximumError([0.5]*3)
data_edge_2 = Canny.Execute(data_32)

可視化

可視化的方法很是簡單 只須要一條指令:

sitk.Show()

但須要先安裝工具ImageJ,不然沒法使用.具體的安裝連接,能夠參考這篇博文:sitk.show()與imageJ結合使用常見的問題

同一張Dicom文件使用sitk.Show()獲得的效果以下圖:
使用ImageJ進行可視化的結果
除此以外,ImageJ還有一個Tool Bar 支持對圖像的進一步處理:

ImageJ Tool Bar
可見,SimpleITK的可視化要比上面介紹的強大不少,不只能夠實現單張圖像的可視化以及圖像處理,還能夠同時對整個序列的圖像進行統一處理.

單張影像的寫入

一樣有兩種方法

  1. sitk.WriteImage()
new_name = "new_MR_2.dcm"
sitk.WriteImage(img,os.path.join(folder_name,new_name))
  1. sitk.ImageFileWriter()
file_writer = sitk.ImageFileWriter()
file_writer.SetFileName(os.path.join(folder_name,new_name))
file_writer.SetImageIO(imageio="GDCMImageIO")
file_writer.Execute(img)

使用這兩種方法進行寫入的時候,會發現,即使什麼也沒有作,但獲得的新Dicom文件要小於原始的Dicom文件.這是由於新的Dicom文件中沒有Private Creator信息(屬於Dicom Tag的內容).固然若是原始Dicom文件中本就沒有這種信息,文件大小是保持相同的.
由於不少時候只是對圖像進行處理,因此再也不深究.


如有錯誤,還請指出,謝謝!

發佈了1 篇原創文章 · 獲贊 2 · 訪問量 98
相關文章
相關標籤/搜索