[转]常见医疗扫描图像处理步骤

一 数据格式

1.1 dicom

DICOM是医学图像中标准文件，这些文件包含了诸多的元数据信息（比如像素尺寸，每个维度的一像素代表真实世界里的长度）。此处以kaggle Data Science Bowl 数据集为例。

data-science-bowl-2017。数据列表如下:

后缀为 .dcm。

每个病人的一次扫描CT(scan)可能有几十到一百多个dcm数据文件(slices)。可以使用 python的dicom包读取，读取示例代码如下:

1. dicom.read_file('/data/lung_competition/stage1/7050f8141e92fa42fd9c471a8b2f50ce/498d16aa2222d76cae1da144ddc59a13.dcm')

其pixl_array包含了真实数据。

1. slices = [dicom.read_file(os.path.join(folder_name,filename)) for filename in os.listdir(folder_name)]
2. slices = np.stack([s.pixel_array for s in slices])

1.2 mhd格式

mhd格式是另外一种数据格式，来源于(LUNA2016)[https://luna16.grand-challenge.org/data/]。每个病人一个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扫描中原点位置的差值，需要做进一步转换。转换步骤参考 SimpleITK图像转换

1.3 查看CT扫描文件软件

一个开源免费的查看软件 mango

二 dicom格式数据处理过程

2.1 处理思路

首先，需要明白的是医学扫描图像(scan)其实是三维图像，使用代码读取之后开源查看不同的切面的切片(slices)，可以从不同轴切割

如下图展示了一个病人CT扫描图中，其中部分切片slices

其次，CT扫描图是包含了所有组织的，如果直接去看，看不到任何有用信息。需要做一些预处理,预处理中一个重要的概念是放射剂量，衡量单位为HU(Hounsfield Unit)，下表是不同放射剂量对应的组织器官

substance	HU
空气	-1000
肺	-500
脂肪	-100到-50
水	0
CSF	15
肾	30
血液	+30到+45
肌肉	+10到+40
灰质	+37到+45
白质	+20到+30
Liver	+40到+60
软组织,constrast	+100到+300
骨头	+700(软质骨)到+3000(皮质骨)

1. Hounsfield Unit = pixel_value * rescale_slope + rescale_intercept

一般情况rescale slope = 1, intercept = -1024。

上表中肺部组织的HU数值为-500,但通常是大于这个值，比如-320、-400。挑选出这些区域，然后做其他变换抽取出肺部像素点。

2.2 先载入必要的包

1. # -*- coding:utf-8 -*-
2. '''
3. this script is used for basic process of lung 2017 in Data Science Bowl
4. '''
6. import glob
7. import os
8. import pandas as pd
9. import SimpleITK as sitk
11. import numpy as np # linear algebra
12. import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
13. import skimage, os
14. from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing
15. from skimage.measure import label,regionprops, perimeter
16. from skimage.morphology import binary_dilation, binary_opening
17. from skimage.filters import roberts, sobel
18. from skimage import measure, feature
19. from skimage.segmentation import clear_border
20. from skimage import data
21. from scipy import ndimage as ndi
22. import matplotlib
23. #matplotlib.use('Agg')
24. import matplotlib.pyplot as plt
25. from mpl_toolkits.mplot3d.art3d import Poly3DCollection
26. import dicom
27. import scipy.misc
28. import numpy as np

如下代码是载入一个扫描面，包含了多个切片(slices)，我们仅简化的将其存储为python列表。数据集中每个目录都是一个扫描面集（一个病人）。有个元数据域丢失，即Z轴方向上的像素尺寸，也即切片的厚度 。所幸，我们可以用其他值推测出来，并加入到元数据中。

1. # Load the scans in given folder path
2. def load_scan(path):
3. slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
4. slices.sort(key = lambda x: int(x.ImagePositionPatient[2]))
5. try:
6. slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
7. except:
8. slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
10. for s in slices:
11. s.SliceThickness = slice_thickness
13. return slices

2.3 灰度值转换为HU单元

首先去除灰度值为 -2000的pixl_array，CT扫描边界之外的灰度值固定为-2000(dicom和mhd都是这个值)。第一步是设定这些值为0，当前对应为空气（值为0）

回到HU单元，乘以rescale比率并加上intercept(存储在扫描面的元数据中)。

1. def get_pixels_hu(slices):
2. image = np.stack([s.pixel_array for s in slices])
3. # Convert to int16 (from sometimes int16),
4. # should be possible as values should always be low enough (<32k)
5. image = image.astype(np.int16)
6. # Set outside-of-scan pixels to 0
7. # The intercept is usually -1024, so air is approximately 0
8. image[image == -2000] = 0
10. # Convert to Hounsfield units (HU)
11. for slice_number in range(len(slices)):
13. intercept = slices[slice_number].RescaleIntercept
14. slope = slices[slice_number].RescaleSlope
16. if slope != 1:
17. image[slice_number] = slope * image[slice_number].astype(np.float64)
18. image[slice_number] = image[slice_number].astype(np.int16)
20. image[slice_number] += np.int16(intercept)
22. return np.array(image, dtype=np.int16)

可以查看病人的扫描HU值分布情况

1. first_patient = load_scan(INPUT_FOLDER + patients[0])
2. first_patient_pixels = get_pixels_hu(first_patient)
3. plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
4. plt.xlabel("Hounsfield Units (HU)")
5. plt.ylabel("Frequency")
6. plt.show()

2.4 重采样

不同扫描面的像素尺寸、粗细粒度是不同的。这不利于我们进行CNN任务，我们可以使用同构采样。

一个扫描面的像素区间可能是[2.5,0.5,0.5],即切片之间的距离为2.5mm。可能另外一个扫描面的范围是[1.5,0.725,0.725]。这可能不利于自动分析。

常见的处理方法是从全数据集中以固定的同构分辨率重新采样，将所有的东西采样为1mmx1mmx1mm像素。

1. def resample(image, scan, new_spacing=[1,1,1]):
2. # Determine current pixel spacing
3. spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))
4. spacing = np.array(list(spacing))
5. resize_factor = spacing / new_spacing
6. new_real_shape = image.shape * resize_factor
7. new_shape = np.round(new_real_shape)
8. real_resize_factor = new_shape / image.shape
9. new_spacing = spacing / real_resize_factor
11. image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
13. return image, new_spacing

现在重新取样病人的像素，将其映射到一个同构分辨率 1mm x1mm x1mm。

1. pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])

使用matplotlib输出肺部扫描的3D图像方法。可能需要一两分钟。

1. def plot_3d(image, threshold=-300):
3. # Position the scan upright,
4. # so the head of the patient would be at the top facing the camera
5. p = image.transpose(2,1,0)
7. verts, faces = measure.marching_cubes(p, threshold)
8. fig = plt.figure(figsize=(10, 10))
9. ax = fig.add_subplot(111, projection='3d')
10. # Fancy indexing: `verts[faces]` to generate a collection of triangles
11. mesh = Poly3DCollection(verts[faces], alpha=0.1)
12. face_color = [0.5, 0.5, 1]
13. mesh.set_facecolor(face_color)
14. ax.add_collection3d(mesh)
15. ax.set_xlim(0, p.shape[0])
16. ax.set_ylim(0, p.shape[1])
17. ax.set_zlim(0, p.shape[2])
18. plt.show()

打印函数有个阈值参数，来打印特定的结构，比如tissue或者骨头。400是一个仅仅打印骨头的阈值(HU对照表)

1. plot_3d(pix_resampled, 400)

2.5 输出一个病人scans中所有切面slices

1. def plot_ct_scan(scan):
2. '''
3. plot a few more images of the slices
4. :param scan:
5. :return:
6. '''
7. f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(50, 50))
8. for i in range(0, scan.shape[0], 5):
9. plots[int(i / 20), int((i % 20) / 5)].axis('off')
10. plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone)

此方法的效果示例如下:

2.6 定义分割出CT切面里面肺部组织的函数

下面的代码使用了pythonde 的图像形态学操作。具体可以参考python高级形态学操作

1. def get_segmented_lungs(im, plot=False):
3. '''
4. This funtion segments the lungs from the given 2D slice.
5. '''
6. if plot == True:
7. f, plots = plt.subplots(8, 1, figsize=(5, 40))
8. '''
9. Step 1: Convert into a binary image.
10. '''
11. binary = im < 604
12. if plot == True:
13. plots[0].axis('off')
14. plots[0].set_title('binary image')
15. plots[0].imshow(binary, cmap=plt.cm.bone)
17. '''
18. Step 2: Remove the blobs connected to the border of the image.
19. '''
20. cleared = clear_border(binary)
21. if plot == True:
22. plots[1].axis('off')
23. plots[1].set_title('after clear border')
24. plots[1].imshow(cleared, cmap=plt.cm.bone)
26. '''
27. Step 3: Label the image.
28. '''
29. label_image = label(cleared)
30. if plot == True:
31. plots[2].axis('off')
32. plots[2].set_title('found all connective graph')
33. plots[2].imshow(label_image, cmap=plt.cm.bone)
34. '''
35. Step 4: Keep the labels with 2 largest areas.
36. '''
37. areas = [r.area for r in regionprops(label_image)]
38. areas.sort()
39. if len(areas) > 2:
40. for region in regionprops(label_image):
41. if region.area < areas[-2]:
42. for coordinates in region.coords:
43. label_image[coordinates[0], coordinates[1]] = 0
44. binary = label_image > 0
45. if plot == True:
46. plots[3].axis('off')
47. plots[3].set_title(' Keep the labels with 2 largest areas')
48. plots[3].imshow(binary, cmap=plt.cm.bone)
49. '''
50. Step 5: Erosion operation with a disk of radius 2. This operation is
51. seperate the lung nodules attached to the blood vessels.
52. '''
53. selem = disk(2)
54. binary = binary_erosion(binary, selem)
55. if plot == True:
56. plots[4].axis('off')
57. plots[4].set_title('seperate the lung nodules attached to the blood vessels')
58. plots[4].imshow(binary, cmap=plt.cm.bone)
59. '''
60. Step 6: Closure operation with a disk of radius 10. This operation is
61. to keep nodules attached to the lung wall.
62. '''
63. selem = disk(10)
64. binary = binary_closing(binary, selem)
65. if plot == True:
66. plots[5].axis('off')
67. plots[5].set_title('keep nodules attached to the lung wall')
68. plots[5].imshow(binary, cmap=plt.cm.bone)
69. '''
70. Step 7: Fill in the small holes inside the binary mask of lungs.
71. '''
72. edges = roberts(binary)
73. binary = ndi.binary_fill_holes(edges)
74. if plot == True:
75. plots[6].axis('off')
76. plots[6].set_title('Fill in the small holes inside the binary mask of lungs')
77. plots[6].imshow(binary, cmap=plt.cm.bone)
78. '''
79. Step 8: Superimpose the binary mask on the input image.
80. '''
81. get_high_vals = binary == 0
82. im[get_high_vals] = 0
83. if plot == True:
84. plots[7].axis('off')
85. plots[7].set_title('Superimpose the binary mask on the input image')
86. plots[7].imshow(im, cmap=plt.cm.bone)
88. return im

此方法每个步骤对图像做不同的处理，依次为二值化、清除边界、连通区域标记、腐蚀操作、闭合运算、孔洞填充、效果如下:

2.7 肺部图像分割

为了减少有问题的空间，我们可以分割肺部图像（有时候是附近的组织）。这包含一些步骤，包括区域增长和形态运算，此时，我们只分析相连组件。

步骤如下：

• 阈值图像（-320HU是个极佳的阈值，但是此方法中不是必要）

• 处理相连的组件，以决定当前患者的空气的标签，以1填充这些二值图像

• 可选：当前扫描的每个轴上的切片，选定最大固态连接的组织（当前患者的肉体和空气），并且其他的为0。以掩码的方式填充肺部结构。

• 只保留最大的气袋（人类躯体内到处都有气袋）

1. def largest_label_volume(im, bg=-1):
2. vals, counts = np.unique(im, return_counts=True)
3. counts = counts[vals != bg]
4. vals = vals[vals != bg]
5. if len(counts) > 0:
6. return vals[np.argmax(counts)]
7. else:
8. return None
9. def segment_lung_mask(image, fill_lung_structures=True):
11. # not actually binary, but 1 and 2.
12. # 0 is treated as background, which we do not want
13. binary_image = np.array(image > -320, dtype=np.int8)+1
14. labels = measure.label(binary_image)
16. # Pick the pixel in the very corner to determine which label is air.
17. # Improvement: Pick multiple background labels from around the patient
18. # More resistant to "trays" on which the patient lays cutting the air
19. # around the person in half
20. background_label = labels[0,0,0]
22. #Fill the air around the person
23. binary_image[background_label == labels] = 2
26. # Method of filling the lung structures (that is superior to something like
27. # morphological closing)
28. if fill_lung_structures:
29. # For every slice we determine the largest solid structure
30. for i, axial_slice in enumerate(binary_image):
31. axial_slice = axial_slice - 1
32. labeling = measure.label(axial_slice)
33. l_max = largest_label_volume(labeling, bg=0)
35. if l_max is not None: #This slice contains some lung
36. binary_image[i][labeling != l_max] = 1
38. binary_image -= 1 #Make the image actual binary
39. binary_image = 1-binary_image # Invert it, lungs are now 1
41. # Remove other air pockets insided body
42. labels = measure.label(binary_image, background=0)
43. l_max = largest_label_volume(labels, bg=0)
44. if l_max is not None: # There are air pockets
45. binary_image[labels != l_max] = 0
47. return binary_image

查看切割效果

1. segmented_lungs = segment_lung_mask(pix_resampled, False)
2. segmented_lungs_fill = segment_lung_mask(pix_resampled, True)
3. plot_3d(segmented_lungs, 0)

我们可以将肺内的结构也包含进来（结节是固体），不仅仅只是肺部内的空气

1. plot_3d(segmented_lungs_fill, 0)

使用mask时，要注意首先进行形态扩充(python的skimage的skimage.morphology)操作（即使用圆形kernel，结节是球体），参考 python形态操作。这会在所有方向（维度）上扩充mask。仅仅肺部的空气+结构将不会包含所有结节，事实上有可能遗漏黏在肺部一侧的结节（这会经常出现，所以建议最好是扩充mask）。

2.8 数据标准化处理

归一化处理

当前的值范围是[-1024,2000]。而任意大于400的值并不是处理肺结节需要考虑，因为它们都是不同反射密度下的骨头。LUNA16竞赛中常用来做归一化处理的阈值集是-1000和400.以下代码

归一化

1. MIN_BOUND = -1000.0
2. MAX_BOUND = 400.0
4. def normalize(image):
5. image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
6. image[image>1] = 1.
7. image[image<0] = 0.
8. return image

0值中心化

简单来说就是所有像素值减去均值。LUNA16竞赛中的均值大约是0.25.

不要对每一张图像做零值中心化（此处像是在kernel中完成的）CT扫描器返回的是校准后的精确HU计量。不会出现普通图像中会出现某些图像低对比度和明亮度的情况

1. PIXEL_MEAN = 0.25
3. def zero_center(image):
4. image = image - PIXEL_MEAN
5. return image

归一化和零值中心化的操作主要是为了后续训练网络，零值中心化是网络收敛的关键。

三 mhd格式数据处理过程

3.1 处理思路

mhd的数据只是格式与dicom不一样，其实质包含的都是病人扫描。处理MHD需要借助SimpleIKT这个包，处理思路详情可以参考Data Science Bowl2017的toturail Data Science Bowl 2017。需要注意的是MHD格式的数据没有HU值，它的值域范围与dicom很不同。

我们以LUNA2016年的数据处理流程为例。参考代码为 LUNA2016数据切割

3.2 载入必要的包

1. import SimpleITK as sitk
2. import numpy as np
3. import csv
4. from glob import glob
5. import pandas as pd
6. file_list=glob(luna_subset_path+"*.mhd")
7. #####################
8. #
9. # Helper function to get rows in data frame associated
10. # with each file
11. def get_filename(case):
12. global file_list
13. for f in file_list:
14. if case in f:
15. return(f)
16. #
17. # The locations of the nodes
18. df_node = pd.read_csv(luna_path+"annotations.csv")
19. df_node["file"] = df_node["seriesuid"].apply(get_filename)
20. df_node = df_node.dropna()
21. #####
22. #
23. # Looping over the image files
24. #
25. fcount = 0
26. for img_file in file_list:
27. print "Getting mask for image file %s" % img_file.replace(luna_subset_path,"")
28. mini_df = df_node[df_node["file"]==img_file] #get all nodules associate with file
29. if len(mini_df)>0: # some files may not have a nodule--skipping those
30. biggest_node = np.argsort(mini_df["diameter_mm"].values)[-1] # just using the biggest node
31. node_x = mini_df["coordX"].values[biggest_node]
32. node_y = mini_df["coordY"].values[biggest_node]
33. node_z = mini_df["coordZ"].values[biggest_node]
34. diam = mini_df["diameter_mm"].values[biggest_node]

3.3 LUNA16的MHD格式数据的值

一直在寻找MHD格式数据的处理方法，对于dicom格式的CT有很多论文根据其HU值域可以轻易地分割肺、骨头、血液等，但是对于MHD没有这样的参考。从LUNA16论坛得到的解释是，LUNA16的MHD数据已经转换为HU值了，不需要再使用slope和intercept来做rescale变换了。此论坛主题下，有人提出MHD格式没有提供pixel spacing(mm) 和 slice thickness(mm) ，而标准文件annotation.csv文件中结节的半径和坐标都是mm单位，最后确认的是MHD格式文件中只保留了体素尺寸以及坐标原点位置，没有保存slice thickness。即，dicom才是原始数据格式。

3.4 坐标体系变换

MHD值的坐标体系是体素，以mm为单位（dicom的值是GV灰度值）。结节的位置是CT scanner坐标轴里面相对原点的mm值，需要将其转换到真实坐标轴位置，可以使用SimpleITK包中的 GetOrigin() ` GetSpacing()`。图像数据是以512x512数组的形式给出的。

坐标变换如下：

相应的代码处理如下:

1. itk_img = sitk.ReadImage(img_file)
2. img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
3. center = np.array([node_x,node_y,node_z]) # nodule center
4. origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm)
5. spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)
6. v_center =np.rint((center-origin)/spacing) # nodule center in voxel space (still x,y,z ordering)

在LUNA16的标注CSV文件中标注了结节中心的X,Y,Z轴坐标，但是实际取值的时候取的是Z轴最后三层的数组(img_array)。

下述代码只提取了包含结节的最后三个slice的数据，代码参考自 LUNA_mask_extraction.py

1. i = 0
2. for i_z in range(int(v_center[2])-1,int(v_center[2])+2):
3. mask = make_mask(center,diam,i_z*spacing[2]+origin[2],width,height,spacing,origin)
4. masks[i] = mask
5. imgs[i] = matrix2int16(img_array[i_z])
6. i+=1
7. np.save(output_path+"images_%d.npy" % (fcount) ,imgs)
8. np.save(output_path+"masks_%d.npy" % (fcount) ,masks)

3.5 查看结节

以下代码用于查看原始CT和结节mask。其实就是用matplotlib打印上一步存储的npy文件。

1. import matplotlib.pyplot as plt
2. imgs = np.load(output_path+'images_0.npy')
3. masks = np.load(output_path+'masks_0.npy')
4. for i in range(len(imgs)):
5. print "image %d" % i
6. fig,ax = plt.subplots(2,2,figsize=[8,8])
7. ax[0,0].imshow(imgs[i],cmap='gray')
8. ax[0,1].imshow(masks[i],cmap='gray')
9. ax[1,0].imshow(imgs[i]*masks[i],cmap='gray')
10. plt.show()
11. raw_input("hit enter to cont : ")

接下来的处理和DICOM格式数据差不多，腐蚀膨胀、连通区域标记等。

参考信息：

灰度值是pixel value经过重重LUT转换得到的用来进行显示的值，而这个转换过程是不可逆的，也就是说，灰度值无法转换为ct值。只能根据窗宽窗位得到一个大概的范围。 pixel value经过modality lut得到Hu，但是怀疑pixelvalue的读取出了问题。dicom文件中存在（0028，0106）（0028，0107）两个tag，分别是最大最小pixel value，可以用来检验你读取的pixel value 矩阵是否正确。

LUT全称look up table，实际上就是一张像素灰度值的映射表，它将实际采样到的像素灰度值经过一定的变换如阈值、反转、二值化、对比度调整、线性变换等，变成了另外一 个与之对应的灰度值，这样可以起到突出图像的有用信息，增强图像的光对比度的作用。