医学图像处理流程

医学图像处理流程
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包来处理。

⽰例代码如下:
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)。