Chenfan Blog

Do cool things that matter.

0%

医学影像 & AI

Posted on Edited on In

嗯，在杭州待了半个月，想了想还是把一个月前的这篇从草稿堆里面恢复出来了。

被老板外派出去参与健培公司的合作项目，正常的话，大概今年8月就能在央视的节目上看到我们的工作啦 ^0^/

这篇日志里面记录的都是网上的开源资料，主要来源于 Kaggle。

参考资料：

Introduction

深度学习的这把火已经烧遍了几乎所有领域，尤其在图像识别方面。

医疗放射图像方面的疾病监测自然也早就有人盯上了。

LUNA16是一个肺部结节检测的数据集…0.0…不知道是关于肺结节的深度学习检测火了才有了这个数据集还是先有了这个数据集然后肺结节检测才火了。

反正在我们这些领域外行来看，好像特别多的机构和公司都在钻这块的研究。其他的，好像有例如乳腺癌检测等等也是个做的比较多的课题。

除了 LUNA16 数据集本身的排名以外，阿里的天池，还有Kaggle都搞了类似的比赛。

不过相较国外的比赛，阿里的最后要求队伍还要上交所有代码，我感觉就有点砸钱套算法的意味了。嗯，你有钱你牛。

Kaggle 上学术分享的氛围很浓，赛方并没有什么要求，而比赛的前几名往往都会主动把代码开源出来供整个社区一起学习。

Pre-Processing

首先是数据的预处理部分。

Reading CT Data

CT 数据扫描出来会保存成 .dcm 格式的数据，每个文件表示一层 CT 扫描的结果。一个病例则是将一组好多层的结果拼在一起的三维影像。

.dcm 格式中包含了像素大小，图片像素与实际物理尺度的关系等等，相当复杂，python 这边可以用 dicom 这个库直接读取 .dcm 文件，处理起来还比较方便。纯 C++ 操作的话，有个叫 simpleitk 的库可以用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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 = '../input/sample_images/'
patients = os.listdir(INPUT_FOLDER)
patients.sort()

上面这段是从 Kaggle 提供的病例数据中把每个病例的目录列在 patients 中。

之后需要做的就是使用 dicom 库把某个病例目录下所有的 .dcm 文件都读进来，然后处理成人眼能够看的灰度图像。

CT 数据用的是一个叫 Hounsfield 的单位，表示身体的某个部分对于 X 光的不透光性（这个好神奇…）。

根据维基百科上的数据：

|Substance|HU|
|-|
|Air|−1000|
|Lung|-700 to −600|
|Fat|−120 to −90|
|Chyle|−30|
|Water|0|
|Urine|-5 to +15|
|Bile|-5 to +15|
|CSF|+15|
|Kidney|+20 to +45|
|Liver|60 ± 6|
|Lymph nodes|+10 to +20|
|Blood|+30 to +45|
|Muscle|+35 to +55|
|White matter|+20 to +30|
|Grey matter|+37 to +45|
|Soft Tissue, Contrast|+100 to +300|
|Bone|+200 (craniofacial bone), +700 (cancellous bone) to +3000 (cortical bone)|

人体中的每个不同组织的透光性不同，因此可以通过这个来直接过滤出需要的部位数据来。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# 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

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()

get_pixels_hu 这个函数返回的直接是一个 numpy 的三维矩阵了，其中每个点的数据就是 CT 拍出来的 HU 值。

CT 扫描出来的各个方向的物理尺度不一定一致，之后需要对原数据进行重新采样。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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

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)

现在我们得到了一个三个维度缩放比例一致的图像，可以输出来看一下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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(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()

plot_3d(pix_resampled, 400)

measure.marching_cubes 这个函数是在三维的图中对二维表面进行分割……或者说就是在找三维凸包？

400 可以从上面的表里面找到是骨骼的 Hounsfield 值，所以这里打印出来的会是一张肋骨和胸骨的图像。

Lung Segmentation

有了数据，后面需要做的是把肺这个部分从整块 CT 图像中分割出来，毕竟后面的处理和操作都是针对肺的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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)

这里肺分割采用的方法是连通域分析。

先把 -320 阈值以上的所有点都筛掉，剩下在这个值以下的就是肺和空气了。

用 measure.label 这个函数做一次连通域分析，包裹在肺外部的是空气，与 [0][0][0] 相连通的部分全部标记成空。

由于 measure.label 函数默认把 0 作为背景值，因此这里用于连通域分析的二值图像的值实际是 1 和 2。

完成第一次标记之后，将数据恢复成 0 和 1 的二值，再做第二次连通域分析。

保留整个图中最大的连通域就是肺了，剩下的部分都是需要被筛掉的空气。

上面这段代码的输出就是题图的结果。

在输出肺的时候用的是 CPU 进行直接渲染，非常消耗内存，显示出来也很慢。

之后可以对这一块初步分割出来的肺部做进一步的处理，让结果更加精准。

例如 Improved Lung Segmentation using Watershed 等等 Kaggle 上的 Kernels 还提供了更多的参考。

To be continued.

Related Posts