介绍
处理这些文件可能是一个挑战，尤其是考虑到它们的异构性。在将这些dicom文件送入你的CNN之前需要进行一些预处理。
幸运的是，我参加了luna16比赛并作为我的大学的计算机辅助诊断课程的一部分，所以我有一定的工作经验去处理这些文件。
在数据送入你的卷积神经网络/其他机器学习方法前，这个教程能为你提供一个有用步骤的全面概述。
我们讲包括什么：
○加载DICOM文件，并添加缺少的元数据。
○将像素值转换为Hounsfield单位（HU），以及这些单位值对应的组织是什么？
○重采样，重新采样到同构分辨率，以消除扫描仪分辨率的差异。
○3D绘图，可视化是非常有用的，看看我们在做什么。
○肺分割
○归一化，很有意义。
○将扫描对齐为零
在我们开始之前，让我们导入一些软件包并确定可用的患者。
%matplotlib inline
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import dicom
import os
import scipy.ndimage
import matplotlib.pyplot as plt
from skimage import measure, morphology
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
# Some constants
INPUT_FOLDER = '/home/zlk/dicom/'
patients = os.listdir(INPUT_FOLDER)
patients.sort()
加载文件
Dicom是医学成像上的标准文件。这是我第一次使用它，但他似乎相当直接。这些文件包含很多元数据（如像素大小，一个像素在现实世界中的每个维度上有多长）.
不同扫描的像素大小还有粗糙度是不同的（例如，两切片之间的距离可能不同），这将会损害CNN方法的性能。我们可以通过同构重采样来处理这个问题，我们稍后会做。
下面是加载扫描的代码，其中包含多个切片，我们将其保存在一个Python列表中。数据集中的每个文件夹都是一次扫描（所以一个病人）。缺少一个元数据字段，即z方向的像素大小，即片厚。幸运的是，我们可以推断这一点，并将其添加到元数据中。
# Load the scans in given folder path
def load_scan(path):
slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
try:
slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
except:
slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
for s in slices:
s.SliceThickness = slice_thickness
return slices
CT扫描的测量单位是Hounsfield单位（HU），它是放射性强度的量度。CT扫描仪经过仔细校准以精确测量HU值。从维基百科：
但是，默认情况下返回的值不在这个单位中。我们来解决这个问题。
一些扫描仪具有圆柱扫描边界，但是输出图像是方形的。落在这些边界之外的像素得到固定值-2000。第一步是将这些值设置为0，这对应于空气。接下来，让我们回到HU单元，通过乘以重新调整的斜率并添加截取（这些截取方便地存储在扫描的元数据中！）。
def get_pixels_hu(slices):
image = np.stack([s.pixel_array for s in slices])
# Convert to int16 (from sometimes int16),
# should be possible as values should always be low enough (<32k)
image = image.astype(np.int16)
# Set outside-of-scan pixels to 0
# The intercept is usually -1024, so air is approximately 0
image[image == -2000] = 0
# Convert to Hounsfield units (HU)
for slice_number in range(len(slices)):
intercept = slices[slice_number].RescaleIntercept
slope = slices[slice_number].RescaleSlope
if slope != 1:
image[slice_number] = slope * image[slice_number].astype(np.float64)
image[slice_number] = image[slice_number].astype(np.int16)
image[slice_number] += np.int16(intercept)
return np.array(image, dtype=np.int16)
我们来看看其中一个病人。
first_patient = load_scan(INPUT_FOLDER + patients[0])
first_patient_pixels = get_pixels_hu(first_patient)
plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
plt.xlabel("Hounsfield Units (HU)")
plt.ylabel("Frequency")
plt.show()
# Show some slice in the middle
plt.imshow(first_patient_pixels[80], cmap=plt.cm.gray)
plt.show()
从维基百科和这个直方图来看，我们可以清楚地看到哪些像素是空气，哪些是组织。在肺部分割中我们将用这个:)
重采样
一个扫描可能有一个像素间距为[2.5, 0.5, 0.5]，这意味着切片之间的距离是2.5毫米。对于不同的扫描，这可能是[1.5, 0.725, 0.725]，这对于自动分析可能是有问题的（例如使用ConvNets）！
处理这个问题的常用方法是将整个数据集重新采样到一定的各向同性分辨率。如果我们选择重新采样一切到1mm1mm1mm像素，我们可以放心的使用3D卷积神经网络，不用学习缩放/切片厚度不变性。这似乎是一个非常简单的步骤，但由于四舍五入，还是有一些边缘情况。而且，这需要相当长的一段时间。下面的代码适合我们（并处理边缘情况）：
def resample(image, scan, new_spacing=[1,1,1]):
# Determine current pixel spacing
spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32)
resize_factor = spacing / new_spacing
new_real_shape = image.shape * resize_factor
new_shape = np.round(new_real_shape)
real_resize_factor = new_shape / image.shape
new_spacing = spacing / real_resize_factor
image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
return image, new_spacing
请注意，当你用这个，保存新的间距！由于四舍五入，这可能会略微偏离所需的间距（上面的脚本选取最好的可能间距舍入）。
让我们重新采样我们的患者的像素为1毫米×1毫米×1毫米的同构分辨率。
pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])
print("Shape before resampling\t", first_patient_pixels.shape)
print("Shape after resampling\t", pix_resampled.shape)
('Shape before resampling\t', (243, 512, 512))
('Shape after resampling\t', (304, 345, 345))
3D绘制扫描
对于可视化来说，能够显示扫描的三维图像是非常有用的。不幸的是，在这个意义上，这个Kaggle docker映像中可用的包是非常有限的，所以我们将使用Marching Cubes算法创建三维对象的近似网格，用matplotlib绘制这个图。相当缓慢和丑陋，但我们能做的最好。
def plot_3d(image, threshold=-300):
# Position the scan upright,
# so the head of the patient would be at the top facing the camera
p = image.transpose(2,1,0)
verts, faces = measure.marching_cubes_classic(p, threshold)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
# Fancy indexing: `verts[faces]` to generate a collection of triangles
mesh = Poly3DCollection(verts[faces], alpha=0.70)
face_color = [0.45, 0.45, 0.75]
mesh.set_facecolor(face_color)
ax.add_collection3d(mesh)
ax.set_xlim(0, p.shape[0])
ax.set_ylim(0, p.shape[1])
ax.set_zlim(0, p.shape[2])
plt.show()
我们的绘图函数需要一个阈值参数，我们可以用它来绘制某些结构，如所有的组织或只有骨头。400是仅显示骨头的良好阈值（参见上面的Hounsfield单位表）。我们开工吧！
plot_3d(pix_resampled,400)
肺分割
为了减少问题的空间，我们可以分割肺部（通常是一些组织）。我和我的学生同事开发的方法相当有效。
这涉及到很多聪明的步骤。它由区域生长和形态操作的一系列应用组成。在这种情况下，我们将只使用连接的组件分析。
步骤：
○图像阈值（-320HU是一个很好的阈值，但这种方法并不重要）
○做连接组件，确定周围空气的标签，在二进制图像中填充1s。
○可选：对于扫描中的每个轴向切片，确定最大的固体连接组件（身体+人周围的空气），并将其他设置为0.这填充了面罩中的肺中的结构。
○只保留最大的气囊（人体内还有其他的空气袋）。
def largest_label_volume(im, bg=-1):
vals, counts = np.unique(im, return_counts=True)
counts = counts[vals != bg]
vals = vals[vals != bg]
if len(counts) > 0:
return vals[np.argmax(counts)]
else:
return None
def segment_lung_mask(image, fill_lung_structures=True):
# not actually binary, but 1 and 2.
# 0 is treated as background, which we do not want
binary_image = np.array(image > -320, dtype=np.int8)+1
labels = measure.label(binary_image)
# Pick the pixel in the very corner to determine which label is air.
# Improvement: Pick multiple background labels from around the patient
# More resistant to "trays" on which the patient lays cutting the air
# around the person in half
background_label = labels[0,0,0]
#Fill the air around the person
binary_image[background_label == labels] = 2
# Method of filling the lung structures (that is superior to something like
# morphological closing)
if fill_lung_structures:
# For every slice we determine the largest solid structure
for i, axial_slice in enumerate(binary_image):
axial_slice = axial_slice - 1
labeling = measure.label(axial_slice)
l_max = largest_label_volume(labeling, bg=0)
if l_max is not None: #This slice contains some lung
binary_image[i][labeling != l_max] = 1
binary_image -= 1 #Make the image actual binary
binary_image = 1-binary_image # Invert it, lungs are now 1
# Remove other air pockets insided body
labels = measure.label(binary_image, background=0)
l_max = largest_label_volume(labels, bg=0)
if l_max is not None: # There are air pockets
binary_image[labels != l_max] = 0
return binary_image
segmented_lungs = segment_lung_mask(pix_resampled, False)
segmented_lungs_fill = segment_lung_mask(pix_resampled, True)
plot_3d(segmented_lungs, 0)
美丽！
但有一件事我们可以解决，可能是一个好主意，包括肺内的结构（因为结节是固体），我们不只是想在肺部空气。
plot_3d(segmented_lungs_fill,0)
那更好。让我们想象一下两者之间的区别。
plot_3d(segmented_lungs_fill-segmented_lungs,0)
很酷，不是吗？无论如何，当你想使用这个面具时，记得首先对它应用一个扩张形态操作（即用一个循环内核）.这样可以在各个方向上扩展蒙版。只有肺部的空气+结构不会包含所有的结节，特别是它会错过那些经常出现在肺的一侧的结节！所以扩大一点面膜
对于某些边缘情况，此分段可能会失败。它依赖于患者外部的空气不与肺中的空气相连的事实。如果患者有气管切开术，情况就不会这样，我不知道这是否存在于数据集中。此外，特别嘈杂的图像（例如由于在下面的图像起搏器）此方法也可能会失败。相反，身体的第二大气囊将被分割。您可以通过检查掩码对应的图像分数来识别这种情况，对于这种情况来说，这将非常小。然后，您可以首先应用形状闭合操作与几毫米大小的内核关闭这些洞，之后它应该工作（或更简单地说，不要使用此图像的掩码）。
归一化
我们的HU值目前从-1024到2000年。超过400的任何东西都不是我们感兴趣的，因为这些都是不同密度的骨头。在LUNA16竞争中常用的一组阈值是-1000和400.下面是一些你可以使用的代码：
MIN_BOUND = -1000.0
MAX_BOUND = 400.0
def normalize(image):
image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
image[image>1] = 1.
image[image<0] = 0.
return image
零中心
作为最后的预处理步骤，建议零中心数据以使您的平均值为0.为此，您只需从所有像素中减去平均像素值即可。
为了确定这个意思，你只需简单地平均整个数据集中的所有图像。如果这听起来像很多工作，我们发现这在LUNA16比赛中大约为0.25
警告：不要将每个图像的平均值置于中心位置（就像这里的一些内核所做的那样）。CT扫描仪被校准以返回准确的HU测量结果。没有像普通图片那样的对比度或亮度较低的图像。
PIXEL_MEAN = 0.25
def zero_center(image):
image = image - PIXEL_MEAN
return image
接下来做什么?
通过这些步骤，您的图像就可以被CNN或其他ML方法使用了:)。您可以离线执行所有这些步骤（一次并保存结果），我建议您这样做，让它在一夜之间运行，因为这可能需要很长时间。
提示：为了节省存储空间，不要预先进行正常化和零居中，但要在网上进行（在培训过程中，加载之后）。如果你不这样做呢，你的形象是int16，均小于float32s容易压缩以及。
