使用Python對(duì)Dicom文件進(jìn)行讀取與寫(xiě)入的實(shí)現(xiàn)
Pydicom
單張影像的讀取
使用 pydicom.dcmread() 函數(shù)進(jìn)行單張影像的讀取,返回一個(gè)pydicom.dataset.FileDataset對(duì)象.
import osimport pydicom# 調(diào)用本地的 dicom file folder_path = r'D:FilesDataMaterials'file_name = 'PA1_0001.dcm'file_path = os.path.join(folder_path,file_name)ds = pydicom.dcmread(file_path)
在一些特殊情況下,比如直接讀取從醫(yī)院拿到的數(shù)據(jù)(未經(jīng)任何處理)時(shí),可能會(huì)發(fā)生以下報(bào)錯(cuò):
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.
可以看到,由于缺失文件元信息頭,無(wú)法直接讀取,只能強(qiáng)行讀取.這種情況可以直接根據(jù)提示,調(diào)整命令為:
ds = pydicom.dcmread(file_path,force=True)
但后續(xù)還會(huì)碰到:
AttributeError: ’Dataset’ object has no attribute ’TransferSyntaxUID’
在網(wǎng)上檢索后發(fā)現(xiàn),可以通過(guò)設(shè)置TransferSyntaxUID來(lái)解決問(wèn)題:
ds.file_meta.TransferSyntaxUID = pydicom.uid.ImplicitVRLittleEndian
這樣就大功告成了(這里實(shí)際上就提前接觸到了下面讀取Dicom Tags的內(nèi)容了)
一些簡(jiǎn)單處理
讀取成功后,我們可以對(duì) Dicom文件 進(jìn)行一些簡(jiǎn)單的處理
讀取并編輯Dicom Tags
可以通過(guò)兩種方法來(lái)讀取Tag的值
使用的Tag的Description
print(ds.PatientID,ds.StudyDate,ds.Modality)
使用 ds.get() 函數(shù). 函數(shù)內(nèi)參數(shù)采用的是Tag ID.幾種簡(jiǎn)單的打開(kāi)Dicom文件的軟件(如RadiAnt DICOM Viewer)都可以直接看到.這里不再贅述.
ds.get(0x00100020) # 這里得到的是PatientID
讀取到相應(yīng)的Tag值后, 也可以將其他的值寫(xiě)入這些Tag.只要最后保存一下就可以了.
借助Numpy與PIL.Image
讀取Dicom文件后,可以借助Numpy以及圖像處理庫(kù)(如PIL.Image)來(lái)進(jìn)行簡(jiǎn)單的處理.
借助Numpy
import numpy as npdata = np.array(ds.pixel_array)
注意這里使用的是 np.array() 而不是 np.asarray(). 因?yàn)榍罢叩母牟⒉粫?huì)帶來(lái)原pixel_array的改變.在轉(zhuǎn)化為ndarray后 可以直接進(jìn)行簡(jiǎn)單的切割和連接,比如截取某一部分和將兩張圖像拼在一起等,之后再寫(xiě)入并保存下來(lái)即可.
借助PIL.Image
from PIL import Imagedata_img = Image.fromarray(ds.pixel_array)data_img_rotated = data_img.rotate(angle=45,resample=Image.BICUBIC,fillcolor=data_img.getpixel((0,0)))
這里展示的是旋轉(zhuǎn), 還有其他功能如resize等.需要注意的是,從Numpy的ndarray轉(zhuǎn)化為Image時(shí),一般會(huì)發(fā)生變化:
print(data.dtype) # int16data_rotated = np.array(data_img_rotated)print(data_img) # int32
只需要指定參數(shù)就可以解決了
data_rotated = np.array(data_img_rotated,dtype = np.int16)
可視化
簡(jiǎn)單的可視化Pydicom沒(méi)有直接的實(shí)現(xiàn)方法,我們可以通過(guò)上面借助Matplotlib以及Image模塊來(lái)實(shí)現(xiàn).但效果有限.
借助 Matplotlib (Pydicom官方文檔中使用的方法)
from matplotlib import pyplotpyplot.imshow(ds.pixel_array,cmap=pyplot.cm.bone)pyplot.show()
效果如圖所示:
但真實(shí)的圖像是:
顯然顏色是有區(qū)別的.導(dǎo)致這種差別的原因是pyplot函數(shù)使用的cm也就是'color map' 是簡(jiǎn)單的'bone' 并不能滿(mǎn)足醫(yī)學(xué)圖像的要求.
借助Image模塊
data_img.show()
一條指令即可,但是效果很差,如圖所示:
綜合來(lái)看,兩種方法都不是很好.
單張影像的寫(xiě)入
經(jīng)過(guò)上面對(duì)Tag值的修改, 對(duì)圖像的切割, 旋轉(zhuǎn)等操作.最后需要重新寫(xiě)入該Dicom文件.
ds.PixelData = data_rotated.tobytes()ds.Rows,ds.Columns = data_rotated.shapenew_name = 'dicom_rotated.dcm'ds.save_as(os.path.join(folder_path,new_name))
SimpleITK
SimpleITK 是從基于C++的ITK遷移到Python的,所以很多方法的使用都跟C++很相似.
import SimpleITK as sitk
單張影像的讀取
有兩種方法:
sitk.ReadImage()這種方法直接返回image對(duì)象,簡(jiǎn)單易懂.但是無(wú)法讀取Tag的值.
img = sitk.ReadImage(file_path)print(type(img)) # <class ’SimpleITK.SimpleITK.Image’>
sitk.ImageFileReader()這種方法比較像C++的操作風(fēng)格,需要先初始化一個(gè)對(duì)象,然后設(shè)置一些參數(shù),最后返回image.相對(duì)更復(fù)雜,但可以操作的點(diǎn)比較多
file_reader = sitk.ImageFileReader()file_reader.SetFileName(file_path) #這里只顯示了必需的,還有很多可以設(shè)置的參數(shù)data = file_reader.Execute()# 使用這種方法讀取Dicom的Tag Valuefor key in file_reader.GetMetaDataKeys(): print(key,file_reader.GetMetaData(key))
以上兩種方法返回的都是三維的對(duì)象,這與Pydicom有很大的不同.
data_np = sitk.GetArrayFromImage(data)print(data_np.shape) # (1, 512, 512) = (Slice index, Rows, Columns)
序列讀取
序列讀取的方法與單張圖像讀取的第二種方法很相似.(暫且只發(fā)現(xiàn)了一種方法讀取序列,如果還有其他方法,請(qǐng)?jiān)谠u(píng)論區(qū)予以補(bǔ)充,感謝!)
series_reader = sitk.ImageSeriesReader()fileNames = series_reader.GetGDCMSeriesFileNames(folder_name)series_reader.SetFileNames(fileNames)images = series_reader.Execute()
同樣,返回的也是三維的對(duì)象.
一些簡(jiǎn)單操作
SimpleITK 包含很多圖像處理如濾波的工具,這里簡(jiǎn)單介紹一個(gè)邊緣檢測(cè)工具和可視化工具
邊緣檢測(cè)
以Canny邊緣檢測(cè)算子為例,與讀取單張圖像類(lèi)似,同樣有兩種方式:
sitk.CannyEdgeDetection()由于濾波的對(duì)象必須是32位圖像或者其他格式, 需要通過(guò) sitk.Cast() 轉(zhuǎn)換. 之后可以再轉(zhuǎn)換回原格式.
data_32 = sitk.Cast(data,sitk.sitkFloat32)data_edge_1 = sitk.CannyEdgeDetection(data_32,5,30,[5]*3,[0.8]*3)
sitk.CannyEdgeDetectionImageFilter()這個(gè)操作相對(duì)麻煩一些
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)
可視化
可視化的方法非常簡(jiǎn)單 只需要一條指令:
sitk.Show()
但需要先安裝工具ImageJ,否則無(wú)法使用.具體的安裝鏈接,可以參考這篇博文:sitk.show()與imageJ結(jié)合使用常見(jiàn)的問(wèn)題
同一張Dicom文件使用sitk.Show()得到的效果如下圖:
除此之外,ImageJ還有一個(gè)Tool Bar 支持對(duì)圖像的進(jìn)一步處理:
可見(jiàn),SimpleITK的可視化要比上面介紹的強(qiáng)大很多,不僅可以實(shí)現(xiàn)單張圖像的可視化以及圖像處理,還可以同時(shí)對(duì)整個(gè)序列的圖像進(jìn)行統(tǒng)一處理.
單張影像的寫(xiě)入
同樣有兩種方法
sitk.WriteImage()
new_name = 'new_MR_2.dcm'sitk.WriteImage(img,os.path.join(folder_name,new_name))
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)
使用這兩種方法進(jìn)行寫(xiě)入的時(shí)候,會(huì)發(fā)現(xiàn),即便什么也沒(méi)有做,但得到的新Dicom文件要小于原始的Dicom文件.這是因?yàn)樾碌腄icom文件中沒(méi)有Private Creator信息(屬于Dicom Tag的內(nèi)容).當(dāng)然如果原始Dicom文件中本就沒(méi)有這種信息,文件大小是保持相同的.因?yàn)楹芏鄷r(shí)候只是對(duì)圖像進(jìn)行處理,所以不再深究.
到此這篇關(guān)于使用Python對(duì)Dicom文件進(jìn)行讀取與寫(xiě)入的實(shí)現(xiàn)的文章就介紹到這了,更多相關(guān)Python Dicom文件進(jìn)行讀取與寫(xiě)入內(nèi)容請(qǐng)搜索好吧啦網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持好吧啦網(wǎng)!
相關(guān)文章:
1. XHTML 1.0:標(biāo)記新的開(kāi)端2. CSS3實(shí)例分享之多重背景的實(shí)現(xiàn)(Multiple backgrounds)3. XML入門(mén)的常見(jiàn)問(wèn)題(四)4. asp(vbscript)中自定義函數(shù)的默認(rèn)參數(shù)實(shí)現(xiàn)代碼5. 詳解瀏覽器的緩存機(jī)制6. php網(wǎng)絡(luò)安全中命令執(zhí)行漏洞的產(chǎn)生及本質(zhì)探究7. XML解析錯(cuò)誤:未組織好 的解決辦法8. 使用Spry輕松將XML數(shù)據(jù)顯示到HTML頁(yè)的方法9. 利用CSS3新特性創(chuàng)建透明邊框三角10. ASP基礎(chǔ)知識(shí)VBScript基本元素講解
