医学图像处理中的数据读写

医学图像处理中的数据读写

常见的医学图像的格式

不管格式如何变化，对于医学图像而言，最终读取到内容中的数据就是图像的强度值信息，就类似自然图像的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 sitkimport sysimport osimport numpy as npimport 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 361Max 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以后，后面的操作都是一样的。

#读取mhamhaimg = sitk.ReadImage(os.path.join(mhadir,"data1.mha"))print(mhaimg.GetSize())print(mhaimg.GetSpacing())#转为numpyrawmhaimg = sitk.GetArrayFromImage(mhaimg)print(rawmhaimg.shape)#保存mhasitk.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']

#读取mhdmhdimg = sitk.ReadImage(os.path.join(mhddir,"data1.mhd"))print(mhdimg.GetSize())print(mhdimg.GetSpacing())#转为numpyrawmhdimg = sitk.GetArrayFromImage(mhdimg)print(rawmhdimg.shape)#保存mhdsitk.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']

#读取niiniiimg = sitk.ReadImage(os.path.join(niidir,"data1.nii.gz"))print(niiimg.GetSize())print(niiimg.GetSpacing())#转为numpyrawniiimg = sitk.GetArrayFromImage(niiimg)print(rawniiimg.shape)#保存niisitk.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']

#读取nrrdnrrdimg = sitk.ReadImage(os.path.join(nrrddir,"data1.nrrd"))print(nrrdimg.GetSize())print(nrrdimg.GetSpacing())#转为numpyrawnrrdimg = sitk.GetArrayFromImage(nrrdimg)print(rawnrrdimg.shape)#保存niisitk.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 spacingmrawimg.SetSpacing((0.58984375, 0.58984375, 0.5))#输出文件，我们将其保存为mhasitk.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)