[Medical Image] 医学图像处理中的数据读写
医学图像处理中的数据读写
常见的医学图像的格式
不管格式如何变化,对于医学图像而言,最终读取到内容中的数据就是图像的强度值信息,就类似自然图像的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)
如果你需要其他的信息,可以查阅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)