医学图像格式预处理

医学图像文件数据存储格式

---

DICOM

DICOM(Digital Imaging and Communications in Medicine)是指医疗数字影像传输协定，是用于医学影像处理、储存、打印、传输的一组通用的标准协定。它包含了文件格式的定义以及网络通信协议。DICOM是以TCP/IP为基础的应用协定，并以TCP/IP联系各个系统。两个能接受DICOM格式的医疗仪器间，可通过DICOM格式的文件，来接收与交换影像及病人资料。

目前，DICOM被广泛应用于放射医疗，心血管成像以及放射诊疗诊断设备（X射线，CT，核磁共振，超声等），并且在眼科和牙科等其它医学领域得到越来越深入广泛的应用。

一个 DICOM 文件包含文件头部和同文件名的*.dcm 图像数据。文件头部的大小取决于它所提供的信息的多少。文件头包含以下信息：病人的 ID，病人的姓名，图像的模态以及其他信息。它定义了帧的数量以及图像的精度。这些信息会被图像浏览器在显示图像时用到。对于一个单词采样，会有很多个 DICOM 文件。

尽管DICOM是MRI采集的标准输出格式，但是，数据分析前往往要把DICOM格式转化为其他分析格式，这主要是因为DICOM数据比较庞大。由于DICOM把每层图像都存储为独立文件，这会导致产生大量较小的数字文件，从而堵塞文件系统，降低分析速度。

dicom 文件读取

使用pydicom包

1.  
   # 安装 Pydicom 模组
2.  
   pip3 install pydicom
3.  
    
4.  
   # 以 whl 文件安装 Pydicom 模组
5.  
   pip3 install pydicom-1.4.1-py2.py3-none-any.whl
6.  
    
7.  
   # 透过 Conda 安装 Pydicom 模组
8.  
   conda install -c conda-forge pydicom

读取数据

1.  
   from pydicom import dcmread
2.  
   from pydicom.data import get_testdata_files
3.  
    
4.  
   # 取得 Pydicom 附带的 DICOM 测试影像路径
5.  
   filename = get_testdata_files('MR_small.dcm')[0]
6.  
    
7.  
   # 读取 DICOM 文件
8.  
   ds = dcmread(filename)
9.  
    
10.  
    # 列出所有元数据（metadata）
11.  
    print(ds)
12.  
    # 这里透过 dcmread 函数读取出来的 ds 是一个数据集（dataset），里面包含了 DICOM 的各种信息以及影像数据。
13.  
     
14.  
    #可以直接取出指定字段的数据，例如取得病人姓名（Patient’s Name）字段的数据：
15.  
    # 输出 Patient Name 数据
16.  
    print(ds.PatientName)

显示数据

1.  
   # Pydicom 为了让读取 DICOM 影像更为方便，
2.  
   #提供了 pixel_array 这个 NumPy 阵列，用户从这里就可以直接取得影像的数据，
3.  
   #不需要处理底层数据格式的问题，配合 matplotlib 即可立即绘制影像数据：
4.  
   import matplotlib.pyplot as plt
5.  
    
6.  
   # 以 matplotlib 绘制影像
7.  
   plt.imshow(ds.pixel_array)
8.  
   plt.show()

编辑影像大小

1.  
   # 缩小影像
2.  
   data_downsampling = ds.pixel_array[::4, ::4]
3.  
    
4.  
   # 将缩小的影像放入原来的 DICOM 数据集
5.  
   ds.PixelData = data_downsampling.tobytes()
6.  
    
7.  
   # 更新影像大小
8.  
   ds.Rows, ds.Columns = data_downsampling.shape
9.  
   # 另存 DICOM 文件
10.  
    ds.save_as("de-identification.dcm")

NRRD

官网：Teem: nrrd

NRRD是一种库和文件格式，旨在支持涉及n维栅格数据的科学可视化和图像处理。 NRRD代表“几乎原始的栅格数据”。除了尺寸的通用性外，NRD还具有相对于类型（8种积分类型，2种浮点类型），书面文件（RAW，ASCII，HEX或GZIP或BZIP2压缩）和Endianness的编码（数据字节顺序是明确的，当类型或编码公开时记录）。除了NRRD格式外，库还可以读取和编写PNG，PPM和PGM图像，以及一些VTK“ bustical_points”数据集。实施了关于栅格数据的大约两打操作，包括量价，切片和裁剪等简单的事情，以及更奇特的东西，例如使用任意分离的内核进行投影，直方图均衡和过滤的重采样（上下和下）。

MITK默认会将医学图像保存为格式为NRRD（Nearly Raw Raster Data）的图像，在这个数据格式中包含：

• 一个单个的数据头文件：为科学可视化和医学图像处理准确地表示N维度的栅格信息。
• 既能分开又能合并的图像文件。

一个NRRD文件的大致格式（带有数据头）如下图所示：

nrrd读取数据

1.  
   import numpy as np
2.  
   import nrrd #pip install pynrrd
3.  
    
4.  
   filename = './test.nrrd'
5.  
   readdata, header = nrrd.read(filename)
6.  
   print(readdata.shape) #保存图片的多维矩阵
7.  
   print(header) #保存图片的相关信息

nrrd转nii格式

1.  
   import os
2.  
   from glob import glob
3.  
   import numpy as np
4.  
   import vtk
5.  
    
6.  
   def readnrrd(filename):
7.  
   """Read image in nrrd format."""
8.  
   reader = vtk.vtkNrrdReader()
9.  
   reader.SetFileName(filename)
10.  
    reader.Update()
11.  
    info = reader.GetInformation()
12.  
    return reader.GetOutput(), info
13.  
     
14.  
    def writenifti(image,filename, info):
15.  
    """Write nifti file."""
16.  
    writer = vtk.vtkNIFTIImageWriter()
17.  
    writer.SetInputData(image)
18.  
    writer.SetFileName(filename)
19.  
    writer.SetInformation(info)
20.  
    writer.Write()
21.  
     
22.  
    if __name__ == '__main__':
23.  
    baseDir = os.path.normpath(r'G:/aurora-nrrd/')
24.  
    files = glob(baseDir '/*label.nrrd')
25.  
    for file in files:
26.  
    m, info = readnrrd(file)
27.  
    writenifti(m, file.replace('label.nrrd', 'label.nii'), info)

使用ITK转换

1.  
   import SimpleITK as sitk
2.  
   data_path=r'./data.nrrd'
3.  
   data,options=nrrd.read(data_path)
4.  
   img = sitk.GetImageFromArray(data)
5.  
   # img = sitk.GetArrayFromImage(img)
6.  
   sitk.WriteImage(img,'data.nii.gz')

NIFTI

大部分医学领域导出dicom格式，但是太复杂了。很多时候，将dicom转换为nifti格式也就是nii格式
一个NIFTI格式主要包含三部分：hdr, ext, img

Nifti 格式最初是为神经影像学发明的。神经影像信息学技术计划（NIFTI）将 NIfTI 格式预设为 ANALYZE7.5 格式的替代品。它最初的应用领域是神经影像，但是也被用在其他领域。这种格式的主要特点就是它包含两个能够将每个体素的索引（i,j,k）和它的空间位置（x,y,z）关联起来的仿射坐标。

DICOM 和 NIfTI 这两种格式的主要区别是：NIfTI 中的图像原始数据被存储成了 3 维图像，而 dicom 一些 2 维的图层。这就使得 NIFTI 更加适合那些应用在 DICOM 上的机器学习的方法，因为它是以 3D 图像建模的。处理一个单独的 NIFTI 文件要比处理成百上千个 dicom 文件更加容易一些。与 DICOM 格式下的好多个文件相比，NIFTI 格式下，每个 3d 图像只有两个文件。

神经成像信息技术创新”将NIFTI格式视为ANALYZE7.5格式的替代品。NIFTI最初是用于神经成像的，但它也适用于一些其他的领域。NIFTI中一个主要的特点在于它包含了两个仿射坐标定义，这两个仿射坐标定义能够将每个立体元素指标（i,j,k）和空间位置（x,y,z）联系起来。 Nibabel是用于读取nifti文件的一个朋友Python库，“oro.nifti”是用于读取nifti数据的一个R工具包。

标准NIfTI图像的扩展名是.nii，包含了头文件及图像资料。由于NIfTI格式和Analyze格式的关系，因此NIfTI格式也可使用独立的图像文件[.img]和头文件[.hdr]。单独的.nii格式文件的优势就是可以用标准的压缩软件[如gzip]，而且一些分析软件包[比如FSL]可以直接读取和写入压缩的.nii文件[扩展名为.nii.gz]。

简而言之，nii格式和.nii.gz格式是一个东西。

hdr/header

这部分数据长度是固定的，当然不同版本可能规定的长度不同，但是同一版本的多个nii文件是相同的。
header里包含的信息有：
--维度，x,y,z，单位是毫米。还有第四个维度，就是时间。这部分储存的主要是四个数字。
--voxel size(体素大小)：毫米单位的x,y,z大小。（也就是spacing）
--数据类型，一般是int16，这个精度不够，最好使用double类型。
--Form和转换矩阵，每一个Form都对应一个转换矩阵。（暂时不知道Form是什么）

Extension

是自己可以随意定义数据的部分，可以自己用。但是通用的软件公司都无法使用这部分。

Image

储存3D或者4D的图像数据

坐标

dicom和nii格式定义了不同的方向，对于nii格式，坐标原点在大脑中某个部位上，方向可以从图上看出。

读取.nii/.nii.gz文件

使用nibabel包读取

1.  
   '''
2.  
   查看和显示nii文件
3.  
   '''
4.  
    
5.  
   import matplotlib
6.  
    
7.  
   matplotlib.use('TkAgg')
8.  
    
9.  
   from matplotlib import pylab as plt
10.  
    import nibabel as nib #使用这个工具
11.  
    from nibabel import nifti1
12.  
    from nibabel.viewers import OrthoSlicer3D
13.  
    # 文件名，nii或nii.gz
14.  
    example_filename = './img/0001.nii.gz'
15.  
     
16.  
    img = nib.load(example_filename)
17.  
    print(img)
18.  
    print(img.header['db_name']) # 输出头信息
19.  
    # shape不一定只有三个参数，打印出来看一下
20.  
    width, height, queue = img.dataobj.shape
21.  
    # 显示3D图像
22.  
    OrthoSlicer3D(img.dataobj).show()
23.  
     
24.  
    num = 1
25.  
    # 按照10的步长，切片，显示2D图像
26.  
    for i in range(0, queue, 10):
27.  
    img_arr = img.dataobj[:, :, i]
28.  
    plt.subplot(5, 4, num)
29.  
    plt.imshow(img_arr, cmap='gray')
30.  
    num = 1
31.  
     
32.  
    plt.show()

使用itk包读取

1.  
   import SimpleITK as sitk
2.  
   import skimage.io as io
3.  
    
4.  
    
5.  
   def read_img(path):
6.  
   img = sitk.ReadImage(path)
7.  
   data = sitk.GetArrayFromImage(img)
8.  
   return data
9.  
    
10.  
     
11.  
    # 显示一个系列图
12.  
    def show_img(data):
13.  
    for i in range(data.shape[0]):
14.  
    io.imshow(data[i, :, :], cmap='gray')
15.  
    print(i)
16.  
    io.show()
17.  
     
18.  
     
19.  
    # 单张显示
20.  
    # def show_img(ori_img):
21.  
    # io.imshow(ori_img[100], cmap='gray')
22.  
    # io.show()
23.  
     
24.  
     
25.  
    if __name__ == "__main__":
26.  
    # window下的文件夹路径
27.  
    path = '..\\md_img\\img\\0001.nii.gz' #路径
28.  
    data = read_img(path)
29.  
    show_img(data)

main函数封装读取

1.  
   import SimpleITK as sitk
2.  
   import numpy as np
3.  
   import os
4.  
    
5.  
   from PIL import Image
6.  
    
7.  
    
8.  
   def read_nii(file_path):
9.  
   ds = sitk.ReadImage(file_path) # 读取nii数据的第一个函数sitk.ReadImage
10.  
    # print('ds: ', ds)
11.  
    data = sitk.GetArrayFromImage(ds) # 把itk.image转为array
12.  
    # print('data: ', data)
13.  
    print('shape_of_data', data.shape)
14.  
    spacing = ds.GetSpacing() # 三维数据的间隔
15.  
    # print('spacing_of_data', spacing)
16.  
    return data
17.  
     
18.  
     
19.  
    # 从十六进制的颜色得到RGB颜色
20.  
    def color(value):
21.  
    digit = list(map(str, range(10))) list("ABCDEF")
22.  
    if isinstance(value, tuple):
23.  
    string = '#'
24.  
    for vi in value:
25.  
    a1 = vi // 16
26.  
    a2 = vi % 16
27.  
    string = digit[a1] digit[a2]
28.  
    return string
29.  
    elif isinstance(value, str):
30.  
    a1 = digit.index(value[1]) * 16 digit.index(value[2])
31.  
    a2 = digit.index(value[3]) * 16 digit.index(value[4])
32.  
    a3 = digit.index(value[5]) * 16 digit.index(value[6])
33.  
    return [a1, a2, a3]
34.  
     
35.  
     
36.  
    def getRGBColor(colorArray):
37.  
    colorMapRGB = []
38.  
    for i in range(len(colorArray)):
39.  
    colorMapRGB.append(color(colorArray[i]))
40.  
    return colorMapRGB
41.  
     
42.  
    # niiDataArray(读取nii文件获得的三维数组)
43.  
    # type(需要得到的图片的类型有1：横断面上下切、2矢状面左右切、3冠状面前后切)
44.  
    # imgId(获得的图片在该面的位置)
45.  
    def getImgFromNiiDataArray(niiDataArray, type, imgId, colorMap):
46.  
    shape = niiDataArray.shape
47.  
    if type == 1:
48.  
    imgData = np.array(niiDataArray[imgId, :, :])
49.  
    elif type == 2:
50.  
    imgData = np.array(niiDataArray[:, imgId, :])
51.  
    elif type == 3:
52.  
    imgData = np.array(niiDataArray[:, :, imgId])
53.  
     
54.  
    imgR = np.zeros(imgData.shape)
55.  
    imgG = np.zeros(imgData.shape)
56.  
    imgB = np.zeros(imgData.shape)
57.  
    imgA = np.zeros(imgData.shape)
58.  
    for i in range(1, 33):
59.  
    imgR[imgData == i] = colorMap[i - 1][0]
60.  
    imgG[imgData == i] = colorMap[i - 1][1]
61.  
    imgB[imgData == i] = colorMap[i - 1][2]
62.  
    imgA[imgData == i] = 255
63.  
     
64.  
    r = Image.fromarray(imgR).convert('L')
65.  
    g = Image.fromarray(imgG).convert('L')
66.  
    b = Image.fromarray(imgB).convert('L')
67.  
    a = Image.fromarray(imgA).convert('L')
68.  
    image = Image.merge('RGBA', (r, g, b, a))
69.  
    return image
70.  
     
71.  
     
72.  
    # 以png格式输出nii文件中三种视图的所有图片
73.  
    # 会在target目录下生成三个文件夹(横断面,矢状面,冠状面)
74.  
    # 里面装了对应得一系列png图片
75.  
    def exportAllImgByPNG(niiDataArray, targetPath):
76.  
    if os.path.exists(os.path.join(targetPath, '横断面')) == False:
77.  
    os.mkdir(os.path.join(targetPath, '横断面'))
78.  
    if os.path.exists(os.path.join(targetPath, '矢状面')) == False:
79.  
    os.mkdir(os.path.join(targetPath, '矢状面'))
80.  
    if os.path.exists(os.path.join(targetPath, '冠状面')) == False:
81.  
    os.mkdir(os.path.join(targetPath, '冠状面'))
82.  
     
83.  
    dataShape = niiDataArray.shape
84.  
    #
85.  
    for i in range(dataShape[0]):
86.  
    imge = getImgFromNiiDataArray(niiDataArray, 1, i, rgbColorMap)
87.  
    imge.save(os.path.join(targetPath, '横断面', '_' str(i) '.png'))
88.  
    for i in range(dataShape[1]):
89.  
    imge = getImgFromNiiDataArray(niiDataArray, 2, i, rgbColorMap)
90.  
    imge.save(os.path.join(targetPath, '矢状面', '_' str(i) '.png'))
91.  
    for i in range(dataShape[2]):
92.  
    imge = getImgFromNiiDataArray(niiDataArray, 3, i, rgbColorMap)
93.  
    imge.save(os.path.join(targetPath, '冠状面', '_' str(i) '.png'))
94.  
     
95.  
     
96.  
    # 33个颜色，分别对应牙槽骨和32颗牙齿的颜色，0为牙槽骨的颜色，1~32为牙齿的颜色
97.  
    colorMap = ['#00AA00', '#F93408', '#F57A34', '#F7951E', '#F6C238', '#FBE92F', '#E5F827', '#B6F313', '#97F922',
98.  
    '#75F72A',
99.  
    '#35F80A', '#23FD23', '#1DF645', '#37F57C', '#03FE8C', '#29FEC4', '#36FBE9', '#2CEBFE', '#29C1FA',
100.  
     '#048BFB',
101.  
     '#1B6CFA', '#153FFE', '#2F2FFD', '#3A10F4', '#610DF4', '#8F0AFD', '#B714F4', '#DE11F3', '#FE30EB',
102.  
     '#FD19BF',
103.  
     '#FB279B', '#FD1C6E', '#FD0532'
104.  
     ]
105.  
      
106.  
      
107.  
      
108.  
      
109.  
     if __name__ =='__main__':
110.  
     rgbColorMap = getRGBColor(colorMap) # 得到各部位颜色
111.  
     # print(rgbColorMap)
112.  
      
113.  
     niifile_path = "./lvsili_label_tooth.nii.gz"
114.  
     # file_path = "./lvsiling_label_alveolar.nii.gz"
115.  
      
116.  
     niiDataArray = read_nii(niifile_path)
117.  
     exportAllImgByPNG(niiDataArray, './teethImg')
118.  
     # imge = getImgFromNiiDataArray(niiDataArray, 1, 100, rgbColorMap)
119.  
     # imge.show()
120.  
     # imge.save('./dog.png')

多个nii.gz文件读取

1.  
   # #导入相关的库
2.  
   import numpy as np
3.  
   import os
4.  
   import nibabel as nib
5.  
   import cv2
6.  
   import pickle
7.  
    
8.  
   # #一些参数的设置
9.  
   # .nii.gz 文件路径
10.  
    nii_training_data_path = "D:/WindowsData/Desktop/test_dataset/dataset/"
11.  
    # 生成的 pickle 文件存放路径
12.  
    pickle_data_saving_path = "D:/WindowsData/Desktop/test_dataset/dataset_pickle/"
13.  
    # 保留非 0 数据比重超过 preserving_ratio 的切片
14.  
    preserving_ratio = 0.25
15.  
     
16.  
    filenames = os.listdir(nii_training_data_path)
17.  
    nii_train_list = []
18.  
    for filename in filenames: # 获取所有.nii.gz的文件名
19.  
    if filename.endswith('nii.gz'):
20.  
    nii_train_list.append(filename)
21.  
    print('.nii.gz文件数为：', len(nii_train_list))
22.  
     
23.  
    nii_train_dataset = []
24.  
    for i, f in enumerate(nii_train_list):
25.  
    data_path = os.path.join(nii_training_data_path, f)
26.  
    data = nib.load(data_path).get_data() # 获取.nii.gz文件中的数据
27.  
    # print('.nii.gz文件的维度为：',data.shape) #宽*高*切片数
28.  
    data = data / np.amax(data) # 所有数据归一化
29.  
    for i in range(data.shape[2]): # 对切片进行循环，选出满足要求的切片
30.  
    img = data[:, :, i] # 每一个切片是一个灰色图像
31.  
    if float(np.count_nonzero(img) / img.size) >= preserving_ratio:
32.  
    img = np.transpose(img, (1, 0)) # 将图片顺时针旋转90度（摆正了）
33.  
    nii_train_dataset.append(img)
34.  
    print('选出符合preserving_ratio的图片有：', len(nii_train_dataset))
35.  
    pickle_train_dataset = np.asarray(nii_train_dataset) # list类型转为数组 切片数*宽*高
36.  
     
37.  
    # 将数据保存为pickle类型，方便读取
38.  
    if not os.path.exists(pickle_data_saving_path):
39.  
    os.makedirs(pickle_data_saving_path)
40.  
    with open(os.path.join(pickle_data_saving_path, 'training.pickle'), 'wb') as f:
41.  
    pickle.dump(pickle_train_dataset, f, protocol=4)
42.  
     
43.  
    print('.nii.gz 转换 .pickle 完成！！！')
44.  
     
45.  
    #######################################################
46.  
    #需要用pickle数据的话，可以用下面方法读取，eg如下
47.  
    f = open(os.path.join(pickle_data_saving_path,'training.pickle'),'rb')
48.  
    train = pickle.load(f) #train 与原来的 pickle_train_dataset 一模一样
49.  
    for i in range(train.shape[0]):
50.  
    cv2.imshow('img',train[i,:,:])
51.  
    cv2.waitKey(100)
52.  
    cv2.destroyAllWindows()

参考资料

1. .nii.gz文件以正确定向灰度图像 - 问答 - Python中文网（.nii.gz文件以正确定向灰度图像）
2. 医学图像——数据读取和预处理 - 灰信网（软件开发博客聚合）（nii,nrrd数据预处理）

这篇好文章是转载于：学新通技术网