医学图像处理流程
Tips: 最近在做医学图像预处理(CT/PET),涉及到了一些盲点和知识点,在这做一些总结。
一、数据格式
-
DICOM
DICOM是医学图像中的标准文件,包含了许多元数据信息,这些信息具体可以分为以下四类:
-
Patient
-
Study
-
Series
-
Image
每一个DICOM Tag都是由两个十六进制数的组合来确定的,分别为Group和Element。如(0010,0010)这个Tag表示的是Patient’s Name,它存储着这张DICOM图像的患者姓名。
每个病人的每个模态的有几十到几百的dcm数据文件(slices)
-
-
mhd
mhd格式是另外一种数据格式。每个病人一个mhd文件和一个同名的raw文件,一个
mhd
通常有几百兆,对应的raw
文件只有1kb。mhd
文件需要借助python的SimpleITK
包来处理。SimpleITK示例代码如下:1. import SimpleITK as sitk 2. itk_img = sitk.ReadImage(img_file) 3. img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering) 4. num_z, height, width = img_array.shape #heightXwidth constitute the transverse plane 5. origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm) 6. spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)
需要注意的是,SimpleITK的
img_array
的数组不是直接的像素值,而是相对于CT扫描中原点位置的差值,需要做进一步转换 -
NIFIT
医学影像早期使用的是DICOM标准,基本上各家厂商都会使用符合DICOM标准的产品,但是这个标准对于数据分析并不方便。在神经影像兴起时就诞生了各种各样的数据存储标准,比如analyze。NIFTI出现的原因是原来一种图像格式是ANALYZE 7.5 format,但是这个图像格式缺少一些信息,比如没有方向信息,病人的左右方位等,如果需要包括额外的信息,就需要一个额外的文件,比如ANALYZE7.5就需要一对<.hdr, .img>文件来保存图像的完整信息。因此,解决这个问题Data Format Working Group (DFWG) 将图像格式完整的定义为NIFTI(Neuroimaging Informatics Technology Initiative)格式。放射学和神经学医生使用软件和使用习惯会有一些不同,所以为了统一格式,DFWG提出NIFTI格式图像。
NIFTI格式的读取软件,也应该可以读取<.hdr, .img>文件对。
NIFTI格式中保存的信息是,一共可以保存7维数据。第1,2,3维度(0维统计使用维度个数)保留给空间维度x,y,z,而第四维度留给时间维度t。第5,6,7维度可供用户其他使用,但是第五维度也有一些预定义的用途,例如存储特定体素分布的参数或存储矢量数据。
二、图像处理
-
dicom数据读取
def readdcm(filepath): #filepath = "./T2" series_id = sitk.ImageSeriesReader.GetGDCMSeriesIDs(filepath) series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(filepath, series_id[0]) series_reader = sitk.ImageSeriesReader() #读# 取数据端口 series_reader.SetFileNames(series_file_names) #读取名称 images = series_reader.Execute()#读取数据 #print(images.GetSpacing()) #sitk.WriteImage(images, "T2_1.nii.gz")#保存为nii return images
-
其他格式
image = sitk.ReadImage(filepath)
三、图像预处理
对于不同的数据类型重采样的方法和目的都不相同。例如在遥感中,重采样是从高分辨率遥感影像中提取出低分辨率影像的过程;在数据挖掘中,重采样是指为了解决训练数据类别不均衡,通过在训练期间通过增加小样本的数量或者减少大样本的数量保持样本类别均衡的算法;在医疗图像中,重采样是指将医疗图像中大小不同的体素归一化到相同的大小。体素是体积元素(Volume Pixel)的简称,一张3D医疗图像可以看成是由若干个体素构成的
不同扫描面的像素尺寸、粗细粒度是不同的。这不利于我们进行CNN任务,我们可以使用同构采样。
一个扫描面的像素区间可能是[2.5,0.5,0.5],即切片之间的距离为2.5mm。可能另外一个扫描面的范围是[1.5,0.725,0.725]。这可能不利于自动分析。
常见的处理方法是从全数据集中以固定的同构分辨率重新采样,将所有的东西采样为1mmx1mmx1mm像素。
对于医疗图像重采样,可以分成三个步骤:首先使用SimpleITK读取数据,获得原始图像的对应的Spacing以及Size,得到图像原始的大小;然后图像原始的大小除以新Spacing得到新Size;最后将新Spacing得到新Size赋值到读取的数据即可。
对于一张大小为128128的彩色图像,其在计算机中可以表示为1281283的矩阵,其中每一个像素点的取值范围为0-255,不同的数值代表不同的亮度。但是对于医疗图像其是由若干个slice组成的,假设每一个slice的大小为512512的单通道的图像,其中每一个像素点表示的是一个体素的取值,其范围可以-1000~2000之间。接下来通过以胰腺分割数据集中PANCREAS_0015.nii.gz为例,对医疗图像中体素这个概念进行讲解。Spacing(0.78125, 0.78125, 1.0)表示的是原始图像体素的大小,也可以将Spacing想象成大小为(0.78125, 0.78125, 1.0)的长方体。而原始图像的Size为 (512, 512, 247),表示的是原始在X轴,Y轴,Z轴中体素的个数。原始图像的size对应的Spacing既可以得到真实3D图像大小(5120.78125,5120.78125,2471 ),在图像重采样只是修改体素的大小,而真实3D图像大小是保持不变的,因此假设我们将Spacing修改成(1.0, 1.0, 2.0)的时候,则修改之后其对应的size应该为((5120.78125)/ 1.0,(5120.78125)/ 1.0,(247*1 ))即(400, 400, 124)。
def transform(image,newSpacing, resamplemethod=sitk.sitkNearestNeighbor):
# 设置一个Filter
resample = sitk.ResampleImageFilter()
# 原来的体素块尺寸
originSize = image.GetSize()
# 原来的体素间的距离
originSpacing = image.GetSpacing()
#print(originSpacing)
# newSize = np.array(newSize, float)
# newSpacing = np.array(newSpacing, float)
# factor = originSpacing / newSpacing
# newSize = originSize * factor
# newSpacing = newSpacing.astype(np.int)
newSize = [
int(np.round(originSize[0] * originSpacing[0] / newSpacing[0])),
int(np.round(originSize[1] * originSpacing[1] / newSpacing[1])),
int(np.round(originSize[2] * originSpacing[2] / newSpacing[2]))
]
# 默认像素值(2)
# resample.SetDefaultPixelValue(image.GetPi);
# 沿着x,y,z,的spacing(3)
# The sampling grid of the output space is specified with the spacing along each dimension and the origin.
resample.SetOutputSpacing(newSpacing)
# 设置original
resample.SetOutputOrigin(image.GetOrigin())
# 设置方向
resample.SetOutputDirection(image.GetDirection())
resample.SetSize(newSize)
# 设置插值方式
resample.SetInterpolator(resamplemethod)
# 设置transform
resample.SetTransform(sitk.Euler3DTransform())
# 默认像素值 resample.SetDefaultPixelValue(image.GetPixelIDValue())
return resample.Execute(image)
if __name__ == '__main__':
for i in range(1, 103):
# Dicom序列所在文件夹路径
strr = str(i).zfill(3)
print("Start_", i, '_converse')
filepath = ''
filepathLabel = ''
image = sitk.ReadImage(filepath) # 读取image文件
# print(image.GetOrigin())
# print(image.GetSpacing())
# print(image.GetSize())
label=sitk.ReadImage(filepathLabel) #读取mask
newImage=transform(image,[1,1,1],sitk.sitkLinear)
newLabel=transform(label,[1,1,1])# 不添加插值方式时,默认为最近邻插值
# mask请用最近邻插值,image用线性插值
name_data = ''# image存储path
name_label='' # label存储path
sitk.WriteImage(newImage, name_data)
sitk.WriteImage(newLabel, name_label)