原文连接https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial/notebook

主要包含以下内容：

• 图像加载
• 像素值转换为HU值
• Resampling（重采样）
• 3D绘图
• 肺分割
• 正则化
• Zero centering

加载图像

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

其中path为文件路径，将path下的文件名拿到后与path拼接，读取结果直接保存在list中。

ImagePositionPatient表示图像的左上角在空间坐标系中的x,y,z坐标,单位是毫米. 如果在检查中,则指该序列中第一张影像左上角的坐标。
SliceLocation表示实际的相对位置，单位为mm。
SliceThicknes表示层厚，单位为mm。

通过ImagePositionPatient的Z方向上的值推测出Z方向上的像素尺寸，也就是切片的厚度。最后只用一个for循环将所有切片的厚度赋值。

在读取文件时遇到了一个错误：

使用dicom读取DICOM文件时出错，原因是文件中有其它格式的数据

像素转化为HU值

实际的CT的像素值并不等于HU值，需要进行转化，公式如下：

HU = pixel_value * RescaleSlope + RescaleIntercept

其中的两个R解释如下：

def get_pixels_hu(slices):
    image = np.stack([s.pixel_array for s in slices])
    image = image.astype(np.int8)
    image[image == -2000] = 0

    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)

先将slices中的数组堆叠起来，然后将value=-2000的值置为0，去除CT扫描边界之外的像素值。
然后按照公式逐个image进行转化。

重采样

扫描得到的切片可能具有不用的间距，导致无法直接送到网络中进行训练，为了处理这个问题，通常采用的办法是将完整数据集重新采样为各个方向具有相同分辨率的图像，此处将其采样为1x1x1mm的像素。

关于重采样，可能会改变图像的像素数量，比如放大图像，通过增加原始图像的像素数量来实现，意味着每个像素的大小是不变的。

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)
    print('spacing:{}'.format(spacing))

    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

根据原始像素大小算出图像的大小，然后按照目标像素大小按比例调整像素数。

比如原始像素大小为0.7，有512个像素，则总长为0.7*512，目标像素大小为2，则新调整后的图像有0.7*512/2个像素，记为Num。那么使用Num/512可以算出要扩展像素的比例。

肺分割

在医学图像中，往往包含很多其它组织，为了减少不必要的影响，可以单独将肺部分割出来。这其中涉及了一些巧妙的步骤，比如区域生长和形态学操作。
步骤如下：
阈值分割（-320是一个好的值）。
对图像进行组件连接，以确定当前患者的CT中空气的标签，并用1进行填充。
可选操作：扫描每个切片，选定最大的连接区域（当前患者的肉体和空气），将其它的设为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

该函数选取最大的连接区域的标签，其中bg表示背景值。当统计的vals中有不等于bg的值，则返回该类别标签的值。

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

先使用阈值将图像分割，然后标记不同连通区域。找到一个边缘的标签（空气，因为在转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值中心化

就是对图像的所有像素值减去均值。我们发现Luna_16比赛的数据均值大约为0.25。

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

PIXEL_MEAN = 0.25

def zero_center(image):
    image = image - PIXEL_MEAN
    return image

为了节省空间，不要预先处理数据，可以在线进行正则化和0值中心化。