在上一节：【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割6（数据预处理） 中，我们已经得到了与mhd图像同seriesUID名称的mask nrrd数据文件了，可以说是一一对应了。

并且，mask的文件，还根据结节被多少人同时标注，区分成了4个文件夹，分别是标注了一、二、三、四次，一共就4个医生参与标注。

再加上官方已经给整理好的肺实质分割的文件，我们就获得了以下这些数据：

1. ct图像数据；
2. 肺实质分割数据；
3. 包含结节位置的mask数据。

一、导言

上述得到的这些，就满足了我们的需求了，都是一一对应的，无论是后续的数据预处理，还是拿过来用于训练，都非常的方便。

但是呢，对于原始的ct数据，他在Z轴上的层厚是不同的，这点可以在dicom文件里面看到，也可以在mhd文件的查询到关于层厚的信息。在这点上，不同的序列，差异是非常大的。表现在一个3维数组的结节上面，在这个维度上就是被压扁，和拉长的样子。

在x和y方向，其实也是存在spacing的差异的，但是这种差异没有像z轴那么夸张的，这里可以选择处理和不处理均可（有些论文进行了处理，有些没有。默认都是512x512大小，resample后会变小）。

至此，本篇的目的就很明确了，是要做下面几件事：

1. 对原始图像进行肺实质提取，将肺区外的部分进行裁剪，或者改为固定像素值；
2. 对图像和结节mask进行resample操作，本篇是zyx均进行resample为1mm。

二、具体实施

怎么做的部分，我们分三部分：

1. 肺实质裁剪
2. image和 nodule mask进行resample操作
3. 获取结节中心点坐标和半径

下面就一一展开

2.1、主函数部分

由于这部分数据量比较多，所以在主函数部分采用了多进程的模式，加快处理速度。需要读进来的数据也就是前面篇章已经处理好的，这里都可以直接使用。

下面就是主函数

import sys
import numpy as np
import scipy.ndimage
from skimage import measure, morphology
import SimpleITK as sitk
from multiprocessing import Pool
import os
import nrrd

###############################################################################
# 将标记的mask，和ct原图，加入左右肺区分割的图像，生成去除noise的，剩下肺区的ct和结节mask
###############################################################################
def main():
    n_consensus = 4
    do_resample = True
    img_dir = './LUNA16/image_combined'
    lung_mask_dir = './LUNA16/seg-lungs-LUNA16'
    nod_mask_dir = os.path.join('./LUNA16/nodule_masks', str(n_consensus))

    save_dir = os.path.join('./LUNA16/preprocessed', str(n_consensus))
    os.makedirs(save_dir, exist_ok=True)

    params_lists = []
    # 多进程处理
    for pid in os.listdir(nod_mask_dir):
        #                         seg-lungs-LUNA16, masks_test/3, seg-lungs-LUNA16, preprocessed_test/3, True
        params_lists.append([pid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample])

    pool = Pool(processes=4)
    pool.map(cropResample_process, params_lists)

    pool.close()
    pool.join()

    pool = Pool(processes=4)
    pool.map(generateBBoxes_label, params_lists)

    pool.close()
    pool.join()

if __name__ == '__main__':
    main()

有两个部分，

• cropResample_process：和名称一样，进行肺实质的crop和resample操作；
• generateBBoxes_label：将处理完毕的结节mask，得到结节中心的坐标和半径。

2.2、肺实质裁剪

这小块的步骤，大概如下：

1. 首先，就是数据读取，这部分的详细介绍，可以参考我之前的这篇文章：【医学影像数据处理】nii、npz、npy、dcm、mhd 的数据格式互转，及多目标分割处理汇总
2. 其次，就是将hu值，转化为0-255的值，也就是函数HU2uint8()，对于这部分，可以参考hu值是如何转为0-255的可视化部分的介绍：【医学影像数据处理】 Dicom 文件格式处理汇总
3. 另外，就是将肺区mask作用到图像上，肺实质外采用pad valud补充
4. 最后，将处理好的image、mask和相关参数存储到本地

代码如下，就该说明的部分都进行注释，相信能轻易看懂。

def load_itk_image(filename):
    """Return img array and [z,y,x]-ordered origin and spacing
    """
    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)

    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))

    return numpyImage, numpyOrigin, numpySpacing

def HU2uint8(image, HU_min=-1200.0, HU_max=600.0, HU_nan=-2000.0):
    """
    Convert HU unit into uint8 values. First bound HU values by predfined min
    and max, and then normalize
    image: 3D numpy array of raw HU values from CT series in [z, y, x] order.
    HU_min: float, min HU value.
    HU_max: float, max HU value.
    HU_nan: float, value for nan in the raw CT image.
    """
    image_new = np.array(image)
    image_new[np.isnan(image_new)] = HU_nan

    # normalize to [0, 1]
    image_new = (image_new - HU_min) / (HU_max - HU_min)
    image_new = np.clip(image_new, 0, 1)
    image_new = (image_new * 255).astype('uint8')

    return image_new

def convex_hull_dilate(binary_mask, dilate_factor=1.5, iterations=10):
    """
    Replace each slice with convex hull of it then dilate. Convex hulls used
    only if it does not increase area by dilate_factor. This applies mainly to
    the inferior slices because inferior surface of lungs is concave.
    binary_mask: 3D binary numpy array with the same shape of the image,
        that only region of interest is True. One side of the lung in this
        specifical case.
    dilate_factor: float, factor of increased area after dilation
    iterations: int, number of iterations for dilation
    return: 3D binary numpy array with the same shape of the image,
        that only region of interest is True. Each binary mask is ROI of one
        side of the lung.
    """
    binary_mask_dilated = np.array(binary_mask)
    for i in range(binary_mask.shape[0]):
        slice_binary = binary_mask[i]

        if np.sum(slice_binary) > 0:
            slice_convex = morphology.convex_hull_image(slice_binary)

            if np.sum(slice_convex) <= dilate_factor * np.sum(slice_binary):
                binary_mask_dilated[i] = slice_convex

    struct = scipy.ndimage.generate_binary_structure(3, 1)
    binary_mask_dilated = scipy.ndimage.binary_dilation(
        binary_mask_dilated, structure=struct, iterations=10)

    return binary_mask_dilated

def apply_lung_mask(image, binary_mask1, binary_mask2, pad_value=170,
                    bone_thred=210, remove_bone=False):
    """
    Apply the binary mask of each lung to the image. Regions out of interest
    are replaced with pad_value.
    image: 3D uint8 numpy array with the same shape of the image.
    binary_mask1: 3D binary numpy array with the same shape of the image,
        that only one side of lung is True.
    binary_mask2: 3D binary numpy array with the same shape of the image,
        that only the other side of lung is True.
    pad_value: int, uint8 value for padding image regions that is not
        interested.
    bone_thred: int, uint8 threahold value for determine parts of image is
        bone.
    return: D uint8 numpy array with the same shape of the image after
        applying the lung mask.
    """
    binary_mask = binary_mask1 + binary_mask2
    binary_mask1_dilated = convex_hull_dilate(binary_mask1)
    binary_mask2_dilated = convex_hull_dilate(binary_mask2)
    binary_mask_dilated = binary_mask1_dilated + binary_mask2_dilated
    binary_mask_extra = binary_mask_dilated ^ binary_mask

    # replace image values outside binary_mask_dilated as pad value
    image_new = image * binary_mask_dilated + pad_value * (1 - binary_mask_dilated).astype('uint8')

    # set bones in extra mask to 170 (ie convert HU > 482 to HU 0;
    # water).
    if remove_bone:
        image_new[image_new * binary_mask_extra > bone_thred] = pad_value

    return image_new

def cropResample_process(params):
    #    seg-lungs-LUNA16, masks_test/3, seg-lungs-LUNA16, preprocessed_test/3, True
    pid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample = params

    print('Preprocessing %s...' % (pid))

    img_org, origin, spacing = load_itk_image(os.path.join(img_dir, '%s.mhd' % (pid)))
    lung_mask, _, _ = load_itk_image(os.path.join(lung_mask_dir, '%s.mhd' % (pid)))
    nodule_mask, _ = nrrd.read(os.path.join(nod_mask_dir, '%s.nrrd' % (pid)))

    # 4-右肺   3-左肺   5-气管
    binary_mask_r = lung_mask == 4
    binary_mask_l = lung_mask == 3
    binary_mask = binary_mask_r + binary_mask_l

    img_org = HU2uint8(img_org)
    img_lungRL = apply_lung_mask(img_org, binary_mask_r, binary_mask_l)

有一个点前面从没有说明过，那就是官方提供的lung mask数组，在这里简要的记录下：

• 数字3，表示左肺
• 数字4，表示右肺
• 数字5，表示气管

还是第一次看到这个按位异或运算符（^），简单的学习了下：

按位异或运算符（^）用于将两个操作数的每个对应位进行逻辑异或操作。如果两个对应位的值相同，则结果为0，否则为1。异或的本质是没有进位的加法。

将dilate膨胀后的binary mask和原始的binary mask求异或运算，对应位的值相同，结果为0，否则为1。那么，得到的结果也就是膨胀出来的那部分，就是bone，这部分在去除bone阶段使用到。

可能会有这样的疑问：为什么不直接image 和lung mask相乘，得到一个分割肺实质后留下来的image呢？反而需要采用凸包优化的方式，多此一举呢？

答：在lung mask里面，肺实质的分割是有误差的。也就是肺实质的分割是沿着肺区边缘的，但是某些结节的位置，恰好在肺区的边界上，且密度很大。那么mask就会呈现一个内凹的一个状态。如果采用上面的方法，这样结节就被抠除了。采用凸包优化，就可以利用稍微扩展肺实质边缘，达到将更多肺区留下来的效果。

但是，对于肺结核等等大病灶的疾病，采用上述取出肺实质的方法就不行。主要是因为肺结核的病种范围比较大，尽管采用了凸包优化，最终还是会切除很大一块肺区位置，这样肺区就不完整了，有些得不偿失。

下面是skimage.morphology.convex_hull_image官方给出的实例，如下：点击直达
作用到我们项目里面，切割后的样子如下：

2.3、resample操作

本篇对resample的操作，在zyx的各个维度上，就雨露均沾，通通调整到1mm的状态，这样得到的一个像素大小，表示的也就是物理大小，不会引起任何一个维度上变形的情况。

代码如下所示：

def resample(image, spacing, new_spacing=[1.0, 1.0, 1.0], order=1):
    """
    Resample image from the original spacing to new_spacing, e.g. 1x1x1
    image: 3D numpy array of raw HU values from CT series in [z, y, x] order.
    spacing: float * 3, raw CT spacing in [z, y, x] order.
    new_spacing: float * 3, new spacing used for resample, typically 1x1x1,
        which means standardizing the raw CT with different spacing all into
        1x1x1 mm.
    order: int, order for resample function scipy.ndimage.interpolation.zoom
    return: 3D binary numpy array with the same shape of the image after,
        resampling. The actual resampling spacing is also returned.
    """
    # shape can only be int, so has to be rounded.
    new_shape = np.round(image.shape * spacing / new_spacing)

    # the actual spacing to resample.
    resample_spacing = spacing * image.shape / new_shape

    resize_factor = new_shape / image.shape

    image_new = scipy.ndimage.zoom(image, resize_factor, mode='nearest', order=order)

    return (image_new, resample_spacing)

if do_resample:
    print('Resampling...')
    img_lungRL, resampled_spacing = resample(img_lungRL, spacing, order=3)
    seg_nod_mask = np.zeros(img_lungRL.shape, dtype=np.uint8)
    for i in range(int(nodule_mask.max())):
        # 一个结节，一个结节的resample
        mask = (nodule_mask == (i + 1)).astype(np.uint8)
        mask, _ = resample(mask, spacing, order=3)
        seg_nod_mask[mask > 0.5] = i + 1

其中在resample函数里面，使用到了scipy.ndimage.zoom操作，直接将原始数据，zoom到新的shape。

scipy.ndimage.zoom(input, zoom, output=None, order=3, mode='constant', cval=0.0, prefilter=True, *, grid_mode=False)[source]

函数中：

• input：The input array
• zoom：The zoom factor along the axes

下面是一段官方案例，展示了zoom前后的变化，可以参考：点击链接直达

from scipy import ndimage, datasets
import matplotlib.pyplot as plt

fig = plt.figure()
ax1 = fig.add_subplot(121)  # left side
ax2 = fig.add_subplot(122)  # right side
ascent = datasets.ascent()
result = ndimage.zoom(ascent, 3.0)
ax1.imshow(ascent, vmin=0, vmax=255)
ax2.imshow(result, vmin=0, vmax=255)
plt.show()

zoom前后的变化，如下所示：

发现这个scipy库还真是好用，后续找时间全面的补充下这块的知识。

2.4、存储到本地

这部分就比较的简单了，主要就是说下数组存储的一些新的：

• npy文件存储一些简单的数组，比如下文的spacing、坐标等等；
• nrrd文件存储多维数组，比如下面的image和mask数组图像，大小是240x320x320大小的；
  以前喜欢用nii作为存储文件，现在发现不太好用，nrrd也可以存储数组，还能存储header头。

下面是代码：

lung_box = get_lung_box(binary_mask, img_lungRL.shape)  # 获取肺区分割的外轮廓

z_min, z_max = lung_box[0]
y_min, y_max = lung_box[1]
x_min, x_max = lung_box[2]

# 裁剪操作，去除肺区外的
img_lungRL = img_lungRL[z_min:z_max, y_min:y_max, x_min:x_max]
if do_resample:
    seg_nod_mask = seg_nod_mask[z_min:z_max, y_min:y_max, x_min:x_max]
else:
    seg_nod_mask = nodule_mask[z_min:z_max, y_min:y_max, x_min:x_max]

np.save(os.path.join(save_dir, '%s_origin.npy' % (pid)), origin)  # origin (3,) 记录三维图像origin坐标信息
if do_resample:
    np.save(os.path.join(save_dir, '%s_spacing.npy' % (pid)), resampled_spacing)  # 记录了resample前后x\y\z三个维度的scale系数
np.save(os.path.join(save_dir, '%s_ebox_origin.npy' % (pid)), np.array((z_min, y_min, x_min)))

nrrd.write(os.path.join(save_dir, '%s_clean.nrrd' % (pid)), img_lungRL)  # 去除掉非肺区后的CT图像
nrrd.write(os.path.join(save_dir, '%s_mask.nrrd' % (pid)), seg_nod_mask)  # 去除掉非肺区后的结节MASK图像

2.5、获取结节中心点坐标和半径

这里获取标记结节的中心点坐标和半径，目的还是为了在裁剪patch等操作时候，能够直接从已经获得的结节里面拿取，直接进行crop操作。

这块的步骤和前面get_lung_box差不多，唯一的区别在于保存下来的是中心点，而不是上面的最大、最小边界坐标。

代码如下：

def generateBBoxes_label(params):
    pid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample = params
    masks, _ = nrrd.read(os.path.join(save_dir, '%s_mask.nrrd' % (pid)))

    bboxes = []
    instance_nums = [num for num in np.unique(masks) if num]
    for i in instance_nums:
        mask = (masks == i).astype(np.uint8)
        zz, yy, xx = np.where(mask)
        d = max(zz.max() - zz.min() + 1, yy.max() - yy.min() + 1, xx.max() - xx.min() + 1)
        bboxes.append(np.array([(zz.max() + zz.min()) / 2., (yy.max() + yy.min()) / 2., (xx.max() + xx.min()) / 2., d]))

    bboxes = np.array(bboxes)
    if not len(bboxes):
        print('%s does not have any nodules!!!' % (pid))

    print('Finished masks to bboxes %s' % (pid))
    np.save(os.path.join(save_dir, '%s_bboxes.npy' % (pid)), bboxes)

三、总结

到这里，本篇内容，结合上一篇的内容，我们对Luna16的数据处理基本上就完成了，也完成了我们最早希望得到的内容：

1. images和mask数组，文件名一一对应；
2. resample操作到1mm；
3. 肺实质外的部分丢弃；

6 和 7 这两个篇章，都是对前面几个章节数据部分的补充，你参考这两篇进行数据处理也行，参考其他的数据处理也行，最终得到的数据形式，只要是一样的就行。