醫學影像處理中的資料讀寫

芒果和小貓發表於2022-03-12

醫學影像處理中的資料讀寫

常見的醫學影像的格式

不管格式如何變化，對於醫學影像而言，最終讀取到內容中的資料就是影像的強度值資訊，就類似自然影像的RGB表示法一樣。這裡叫做強度值，因為不同的醫學影像例如CT、MRI他們可區分的訊號的範圍和一般自然影像不一樣。一般自然影像256個級別就夠用了，而醫學影像可能會有2000個級別。對於我而言，比較熟悉的有幾種常見的資料字尾。

  • dicom：常見於臨床上醫學裝置採集的資訊，裡面有許多欄位，我們只關心影像相關的資訊即可。
  • mha：是一種體資料的儲存格式，由一個描述資料的頭和資料組成，一般就針對醫學影像使用。
  • mhd：和mha類似，只不過是把頭和影像資料分開儲存了。
  • nii：NIFTI格式，常用於神經影像，有時候也會見到hdr img這樣的檔案，本質上也是這種格式，只不過是把頭和資料分開了。
  • nrrd：也是一樣包含了頭和影像資料。
  • raw：一般只有影像資訊，其他資訊需要去找對應的解析或者頭。
  • img：一般只有影像資訊，其他資訊需要去找對應的解析或者頭。
    本文沒有過多的介紹這些格式，因為大家主要的目的是做醫學影像處理，其實也不關心他們的具體實現，可以認為他們都是一樣的，都儲存了醫學影像（也可能是三維影像）資訊，我們只需要呼叫現成的庫，將他們解析成類似三維陣列的格式，就可以供我們之後的使用了。

在Python中進行讀寫

Python提供了豐富的庫來輔助我們工作，大多數情況下只需要幾行程式碼就可以搞定讀寫。千萬別想不去去用C++之類的語言，只有在需要工程化的時候再考慮去用C++。

DICOM的讀寫

我從這個網站 https://www.dicomlibrary.com/ 下載了DICOM Samples DICOM的資料作為Example。

import SimpleITK as sitk
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

#本文從上述網站下載的dicom是一個序列有300+切片資料
dicomdir = "Dataset/dicom/data1/series-00000"
reader = sitk.ImageSeriesReader()

dicom_names = reader.GetGDCMSeriesFileNames(dicomdir)
reader.SetFileNames(dicom_names)

image = reader.Execute()

size = image.GetSize()
print("Image size:", size[0], size[1], size[2])

rawImg = sitk.GetArrayFromImage(image)
print("Max value:",np.max(rawImg),"Min value:",np.min(rawImg),rawImg.shape)

Image size: 512 512 361
Max value: 2948 Min value: -1000 (361, 512, 512)

#以上程式碼成功讀取了dicom序列檔案，下面視覺化幾個切片看看
imgslicer = []
for i in range(4):
    #擴充維度，為了視覺化
    s = np.expand_dims(rawImg[i*70,:,:],2)
    #增加通道
    s = np.repeat(s,3,2)
    #這裡進行歸一化，因為dicom這裡有3000個level所以歸一化肯定會丟失精度
    #一般的做法是給定一個區間，再歸一化，其他區域就不管了，這個叫開窗
    print(s.shape)
    s = s.astype(np.float32)
    s = (s - np.min(s))/(np.max(s) - np.min(s))
    imgslicer.append(s)

# rawImg = rawImg.transpose((1,2,0))
# rawImg = np.repeat(rawImg,3,2)
# print("Max value:",np.max(rawImg),"Min value:",np.min(rawImg),rawImg.shape)
for i in range(4):
    plt.subplot(1,4,i+1)
    plt.imshow(imgslicer[i])

(512, 512, 3)
(512, 512, 3)
(512, 512, 3)
(512, 512, 3)

img

如果你需要其他的資訊，可以查閱SimpleITK的文件，另外寫入dicom不在這裡說明了。

mha檔案的讀寫

mha檔案格式更加簡單和緊湊，也是個人用的比較多的一種格式，現在我們將上面使用的dicom檔案轉成mha檔案，來作為我們的Example。注意，有些重複的程式碼這裡就不寫了，包括import庫等，如果後面需要新的庫，會重新import。

mhadir = os.path.join("Dataset","mha")
if not os.path.exists(mhadir):
    os.mkdir(mhadir)
sitk.WriteImage(image,os.path.join(mhadir,"data1.mha"))

在得到了mha檔案以後我們就可以讀取mha檔案了,這裡就不進行視覺化，因為轉成numpy以後，後面的操作都是一樣的。

#讀取mha
mhaimg = sitk.ReadImage(os.path.join(mhadir,"data1.mha"))
print(mhaimg.GetSize())
print(mhaimg.GetSpacing())

#轉為numpy
rawmhaimg = sitk.GetArrayFromImage(mhaimg)
print(rawmhaimg.shape)

#儲存mha
sitk.WriteImage(mhaimg,os.path.join(mhadir,"data1_cpy.mha"))

(512, 512, 361)
(0.58984375, 0.58984375, 0.5)
(361, 512, 512)

mhd檔案的讀寫

mhd檔案類似於mha檔案，只不過是將檔案分開進行了儲存，我們繼續將之前的檔案轉為mhd作為Example

mhddir = os.path.join("Dataset","mhd")
if not os.path.exists(mhddir):
    os.mkdir(mhddir)
sitk.WriteImage(image,os.path.join(mhddir,"data1.mhd"))
print(os.listdir(mhddir))

['data1.mhd', 'data1.raw']

#讀取mhd
mhdimg = sitk.ReadImage(os.path.join(mhddir,"data1.mhd"))
print(mhdimg.GetSize())
print(mhdimg.GetSpacing())

#轉為numpy
rawmhdimg = sitk.GetArrayFromImage(mhdimg)
print(rawmhdimg.shape)

#儲存mhd
sitk.WriteImage(mhaimg,os.path.join(mhddir,"data1_cpy.mhd"))

(512, 512, 361)
(0.58984375, 0.58984375, 0.5)
(361, 512, 512)

nii檔案的讀寫

繼續使用dicom讀取的資料進行讀寫說明

niidir = os.path.join("Dataset","nii")
if not os.path.exists(niidir):
    os.mkdir(niidir)
sitk.WriteImage(image,os.path.join(niidir,"data1.nii.gz"))
print(os.listdir(niidir))

['data1.nii.gz']

#讀取nii
niiimg = sitk.ReadImage(os.path.join(niidir,"data1.nii.gz"))
print(niiimg.GetSize())
print(niiimg.GetSpacing())

#轉為numpy
rawniiimg = sitk.GetArrayFromImage(niiimg)
print(rawniiimg.shape)

#儲存nii
sitk.WriteImage(niiimg,os.path.join(niidir,"data1_cpy.nii.gz"))

(512, 512, 361)
(0.58984375, 0.58984375, 0.5)
(361, 512, 512)

nrrd檔案讀寫

nrrddir = os.path.join("Dataset","nrrd")
if not os.path.exists(nrrddir):
    os.mkdir(nrrddir)
sitk.WriteImage(image,os.path.join(nrrddir,"data1.nrrd"))
print(os.listdir(nrrddir))

['data1.nrrd']

#讀取nrrd
nrrdimg = sitk.ReadImage(os.path.join(nrrddir,"data1.nrrd"))
print(nrrdimg.GetSize())
print(nrrdimg.GetSpacing())

#轉為numpy
rawnrrdimg = sitk.GetArrayFromImage(nrrdimg)
print(rawnrrdimg.shape)

#儲存nii
sitk.WriteImage(nrrdimg,os.path.join(nrrddir,"data1_cpy.nrrd"))

(512, 512, 361)
(0.58984375, 0.58984375, 0.5)
(361, 512, 512)

raw/img檔案的讀取

關於raw檔案的讀取可能會相對來說要麻煩一些，因為raw格式只存資料，至於如何解析，需要額外的說明，這個說明可能就是像mhd那樣有一個檔案來說明，或者有些是在網站或者文件裡直接給出的，需要自己簡單處理一下，之間就碰到過，拿到的資料是raw/img的字尾，只有純資料資訊，而關於尺寸、資料型別等資訊都是這個資料的網站上給出的，因此需要自己簡單處理一下。本文可以利用mhd格式中那個raw檔案來說明一下。

#我們來看一下之前那個資料的型別是什麼，
#隨便找一個numpy陣列看一下
print(rawmhaimg.dtype)

int16

#利用np從檔案讀取檔案
rawpath = os.path.join(mhddir,"data1.raw")
mdata = np.fromfile(rawpath,dtype=np.int16)
#reshape資料，這裡要注意讀入numpy的時候，對應是(z,y,x)
mdata = mdata.reshape((361, 512, 512))
mrawimg = sitk.GetImageFromArray(mdata)
#設定pixel spacing
mrawimg.SetSpacing((0.58984375, 0.58984375, 0.5))
#輸出檔案，我們將其儲存為mha
sitk.WriteImage(mrawimg,os.path.join(mhadir,"data1_raw2mha.mha"))
#讀取一下看看
rawtest = sitk.ReadImage(os.path.join(mhadir,"data1_raw2mha.mha"))
print(rawtest.GetSize())
print(rawtest.GetSpacing())

(512, 512, 361)
(0.58984375, 0.58984375, 0.5)

相關文章