Candidate Generation and LUNA16 preprocessing

在這個kernel中，我們將討論有助於更好地理解問題陳述和數據可視化的方法。 我還將提供有用的資源和信息的鏈接。

此腳本是用Python編寫的。 我建議人們在桌面上安裝anaconda，因為here提到了它的優點。 本教程中用於讀取，處理和可視化數據的庫是matplotlib，numpy，skimage和pydicom.。

圖像大小（z，512,512），其中z是CT掃描中的切片數量，取決於掃描儀的分辨率。 由於計算能力的限制，這樣的大圖像不能直接送到卷積網絡中。 因此，我們將不得不找到更可能患有癌症的地區。 我們將通過首先分割肺部，然后去除低強度區域來縮小我們的搜索空間。

在本教程中，我們將首先讀入數據集並對其進行可視化。 之后，我們將分割肺部結構，然后使用圖像處理方法在CT掃描中找到感興趣的區域（可能的癌症區域）。 接下來我將介紹如何預處理LUNA16數據集，用於UNet等訓練架構的分割和候選分類。

肺結構的分割是非常具有挑戰性的問題，因為在肺部區域不存在同質性，在肺部結構，不同掃描儀和掃描協議中密度相似。 分割的肺可以進一步用於發現肺部結節候選者和感興趣的區域，這可能有助於更好地分類CT掃描。 發現肺結節區域是一個非常困難的問題，因為附着在血管上的結節或者存在於肺區域的邊界處。 肺結節候選者可以通過切割周圍的三維體素並通過可以在LUNA16數據集上訓練的3D CNNs進一步用於分類。 LUNA 16數據集具有每個CT掃描中的結節的位置，因此對於訓練分類器將是有用的。

 

Reading a CT Scan

輸入文件夾有三件事，一個是sample_images文件夾，它具有CT掃描樣本。 stage1_labels.csv包含階段1訓練集圖像的癌症基礎事實，stage1_sample_submission.csv顯示階段1的提交格式。

In [1]:

# This Python 3 environment comes with many helpful analytics libraries installed
# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python
# For example, here's several helpful packages to load in 

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import skimage, os
from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing
from skimage.measure import label,regionprops, perimeter
from skimage.morphology import binary_dilation, binary_opening
from skimage.filters import roberts, sobel
from skimage import measure, feature
from skimage.segmentation import clear_border
from skimage import data
from scipy import ndimage as ndi
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import dicom
import scipy.misc
import numpy as np

# Input data files are available in the "../input/" directory.
# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory

from subprocess import check_output
print(check_output(["ls", "../input"]).decode("utf8"))
#上面兩句話是無法在windows下執行的（沒有ls的binay文件）
import os

cwd = os.getcwd()
file_path=cwd+'/input'
print(os.listdir(file_path))

結果：

['sample_images', 'stage1_labels.csv', 'stage1_sample_submission.csv']

 

每個三維CT掃描由多個切片組成，其數目取決於掃描儀的分辨率，每個切片都有一個與之關聯的實例號，它告訴從頂部切片的索引。 CT掃描的所有dicom文件都在CT掃描名稱的一個文件夾內。 現在我們將讀取掃描的所有dicom切片，然后根據它們的實例編號堆疊它們以獲得3D肺部CT掃描圖像。

sample_images文件夾大約有20個文件夾，每個文件夾對應一個CT掃描。 在文件夾里面有很多dicom文件。

In [2]:

image_path=file_path+'/sample_images'
print(os.listdir(image_path))

結果：

['00cba091fa4ad62cc3200a657aeb957e', '0a099f2549429d29b32f349e95fb2244', '0a0c32c9e08cc2ea76a71649de56be6d', '0a38e7597ca26f9374f8ea2770ba870d', '0acbebb8d463b4b9ca88cf38431aac69', '0b20184e0cd497028bdd155d9fb42dc9', '0bd0e3056cbf23a1cb7f0f0b18446068', '0c0de3749d4fe175b7a5098b060982a1', '0c37613214faddf8701ca41e6d43f56e', '0c59313f52304e25d5a7dcf9877633b1', '0c60f4b87afcb3e2dfa65abbbf3ef2f9', '0c98fcb55e3f36d0c2b6507f62f4c5f1', '0c9d8314f9c69840e25febabb1229fa4', '0ca943d821204ceb089510f836a367fd', '0d06d764d3c07572074d468b4cff954f', '0d19f1c627df49eb223771c28548350e', '0d2fcf787026fece4e57be167d079383', '0d941a3ad6c889ac451caf89c46cb92a', '0ddeb08e9c97227853422bd71a2a695e', '0de72529c30fe642bc60dcb75c87f6bd']

 

每個CT掃描由多個以DICOM格式提供的2D切片組成。 首先，我將讀取CT掃描的隨機dicom文件。 讀取圖像文件后，我們將用0更新-2000的亮度值，因為它們是落在掃描儀邊界之外的像素。

# Any results you write to the current directory are saved as output.
def loadFileInformation(filename):
    information = {}
    ds = dicom.read_file(filename)
    information['PatientID'] = ds.PatientID
    information['PatientName'] = ds.PatientName
    information['PatientBirthDate'] = ds.PatientBirthDate
    return information

random_dcm_path=image_path+'/00cba091fa4ad62cc3200a657aeb957e/38c4ff5d36b5a6b6dc025435d62a143d.dcm'
# print(loadFileInformation(random_dcm_path))
lung = dicom.read_file(random_dcm_path)

slice = lung.pixel_array
slice[slice == -2000] = 0
plt.imshow(slice, cmap=plt.cm.gray)
plt.show()

結果：

def read_ct_scan(folder_name):
    # Read the slices from the dicom file
    slices = [dicom.read_file(folder_name + '/'+filename) for filename in os.listdir(folder_name)]

    # Sort the dicom slices in their respective order
    slices.sort(key=lambda x: int(x.InstanceNumber))

    # Get the pixel values for all the slices
    slices = np.stack([s.pixel_array for s in slices])
    slices[slices == -2000] = 0
    return slices
某個病人CT_path=image_path+'/00cba091fa4ad62cc3200a657aeb957e'
ct_scan = read_ct_scan(某個病人CT_path)

為了可視化切片，我們將必須繪制它們。 matplotlib用於繪制切片。 plot_ct_scan函數將三維CT掃描圖像數組作為輸入，並繪制等間距切片。 CT掃描是灰度圖像，即每個像素的值是單個樣本，這意味着它僅攜帶強度信息。

def plot_ct_scan(scan):
    f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(25, 25))
    for i in range(0, scan.shape[0], 5):
        plots[int(i / 20), int((i % 20) / 5)].axis('off')
        plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone)
plot_ct_scan(ct_scan)
plt.show()

結果：

 

Segmentation of Lungs

在閱讀CT掃描之后，預處理的第一步就是肺結構的分割，因為顯然感興趣的區域位於肺內。 可以看出，肺部是CT掃描中較暗的區域。 肺內的明亮區域是血管或空氣。 在所有地方使用604（-400HU）的閾值，因為在實驗中發現它工作得很好。 我們從CT掃描圖像的每個切片中分割出肺部結構，盡量不要放過附着在肺壁上的可能區域。 有一些可能附着在肺壁上的結節。

我將首先解釋一個常見的方法，使用簡單的圖像處理和形態學操作來分割肺部，然后將給文獻的良好鏈接提供參考和總結。

def get_segmented_lungs(im, plot=False):
    '''
    This funtion segments the lungs from the given 2D slice.
    '''
    if plot == True:
        f, plots = plt.subplots(8, 1, figsize=(5, 40))
    '''
    Step 1: Convert into a binary image.
    '''
    binary = im < 604
    if plot == True:
        plots[0].axis('off')
        plots[0].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 2: Remove the blobs connected to the border of the image.
    '''
    cleared = clear_border(binary)
    if plot == True:
        plots[1].axis('off')
        plots[1].imshow(cleared, cmap=plt.cm.bone)
    '''
    Step 3: Label the image.
    '''
    label_image = label(cleared)
    if plot == True:
        plots[2].axis('off')
        plots[2].imshow(label_image, cmap=plt.cm.bone)
    '''
    Step 4: Keep the labels with 2 largest areas.
    '''
    areas = [r.area for r in regionprops(label_image)]
    areas.sort()
    if len(areas) > 2:
        for region in regionprops(label_image):
            if region.area < areas[-2]:
                for coordinates in region.coords:
                    label_image[coordinates[0], coordinates[1]] = 0
    binary = label_image > 0
    if plot == True:
        plots[3].axis('off')
        plots[3].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 5: Erosion operation with a disk of radius 2. This operation is
    seperate the lung nodules attached to the blood vessels.
    '''
    selem = disk(2)
    binary = binary_erosion(binary, selem)
    if plot == True:
        plots[4].axis('off')
        plots[4].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 6: Closure operation with a disk of radius 10. This operation is
    to keep nodules attached to the lung wall.
    '''
    selem = disk(10)
    binary = binary_closing(binary, selem)
    if plot == True:
        plots[5].axis('off')
        plots[5].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 7: Fill in the small holes inside the binary mask of lungs.
    '''
    edges = roberts(binary)
    binary = ndi.binary_fill_holes(edges)
    if plot == True:
        plots[6].axis('off')
        plots[6].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 8: Superimpose the binary mask on the input image.
    '''
    get_high_vals = binary == 0
    im[get_high_vals] = 0
    if plot == True:
        plots[7].axis('off')
        plots[7].imshow(im, cmap=plt.cm.bone)

    return im

get_segmented_lungs函數分割CT掃描的二維切片。 為了更好的可視化和理解代碼和應用操作，我已經輸出切片。

get_segmented_lungs(ct_scan[71], True)
plt.show()

結果：

 

 

現在，我將逐片分割整個CT掃描並顯示CT掃描的一些切片。

def segment_lung_from_ct_scan(ct_scan):
    return np.asarray([get_segmented_lungs(slice) for slice in ct_scan])

 

segmented_ct_scan = segment_lung_from_ct_scan(ct_scan)
plot_ct_scan(segmented_ct_scan)
plt.show()

結果：

 

Nodule Candidate/Region of Interest Generation

在從CT掃描圖像中分割出肺部結構之后，我們的任務是找到具有結節的候選區域，因為搜索空間非常大。 而且由於計算的局限性，整幅圖像不能直接用3D CNNs進行分類，需要找到可能的癌症區域並對其進行分類。 在實驗中發現，所有的興趣區域的強度> 604（-400 HU）。 所以，我們使用這個閾值來過濾較暗的區域。 這大大減少了候選人數量，保留了召回率高的重要區域。 然后我們將所有候選點分類以減少誤報。

 

segmented_ct_scan[segmented_ct_scan < 604] = 0
plot_ct_scan(segmented_ct_scan)
plt.show()

 

結果：

 

過濾后，由於血管仍然有很多噪音。 因此，我們進一步刪除兩個最大的連接組件（ remove the two largest connected component）。

selem = ball(2)
binary = binary_closing(segmented_ct_scan, selem)

label_scan = label(binary)

areas = [r.area for r in regionprops(label_scan)]
areas.sort()

for r in regionprops(label_scan):
    max_x, max_y, max_z = 0, 0, 0
    min_x, min_y, min_z = 1000, 1000, 1000

    for c in r.coords:
        max_z = max(c[0], max_z)
        max_y = max(c[1], max_y)
        max_x = max(c[2], max_x)

        min_z = min(c[0], min_z)
        min_y = min(c[1], min_y)
        min_x = min(c[2], min_x)
    if (min_z == max_z or min_y == max_y or min_x == max_x or r.area > areas[-3]):
        for c in r.coords:
            segmented_ct_scan[c[0], c[1], c[2]] = 0
    else:
        index = (max((max_x - min_x), (max_y - min_y), (max_z - min_z))) / (min((max_x - min_x), (max_y - min_y), (max_z - min_z)))

plot_3d函數繪制CT掃描的3D numpy數組。

In [14]:

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)
    p = p[:, :, ::-1]

    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.1)
    face_color = [0.5, 0.5, 1]
    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])

In [15]:

plot_3d(segmented_ct_scan, 604)
plt.show()

結果：

 

還有更多的方法可以去嘗試獲取感興趣的區域。

　　1. 因為我們知道結節是球形的，所以左側血管可以進一步使用形狀特性進行過濾。

　　2. 可以進一步使用超像素（應該指的是加強分倍率）分割，並且可以在分割區域上應用形狀屬性。
　　3. 像UNet這樣的CNN體系結構也可以用來生成候選感興趣區域。

使用圖像處理方法也有很多關於肺分割和結節候選生成的優秀論文。 我列出一些方法的小結：

1.  Automatic segmentation of lung nodules with growing neural gas and support vector machine-所提出的方法包括采集肺的計算機斷層掃描圖像（CT），通過提取胸腔，提取肺和重建實質的原始形狀的技術減少感興趣的體積。 之后，使用生長的神經氣體（GNG）來限制比肺實質更密集的結構（結節，血管，支氣管等）。 下一個階段是類似肺結節的結構與其他結構如血管和支氣管的分離。 最后，通過形狀和紋理測量以及支持向量機將結構分類為結核或非結核。
2. Automated Segmentation of Lung Regions using Morphological Operators in CT scan-get_segmented_lungs方法是這篇文章的一個小修改。 本文采用了604的門檻值。
3. Pre-processing methods for nodule detection in lung CT- 本文采用點增強濾波器應用於體素數據的三維矩陣。 該3D濾波器試圖確定每個體素的局部幾何特性，計算Hessian矩陣的特征值並評估特意構建的“似然”函數，以區分線性，平面和球形物體的局部形態，將其建模為具有3D高斯部分 （Q. Li，S. Sone和K. Doi [6]）。 通過將這種3D濾波器應用於人造圖像，即使在將其疊加到非高斯區域的情況下，我們也驗證了高斯樣區域檢測的效率。
4. Lung Nodule Detection using a Neural Classifier:本文討論了用於結節候選選擇的點增強濾波器和用於假陽性發現減少的神經分類器。 該性能被評估為全自動計算機化的方法，用於在篩查CT中識別可能在視覺解釋期間遺漏的肺癌的肺結節的檢測。

接下來我將討論關於LUNA16數據集的預處理，並用它來訓練UNet模型。

 

UNET for Candidate Point Generation

目前深度學習方法在醫學影像分割問題上取得了較好的效果。 一個非常有名的建築是UNET，在我們的例子中可以用於結節候選點生成。 這些網絡的培訓使用注釋數據集完成。 上述用於候選點生成的圖像處理方法不需要任何訓練數據。 我們使用LUNA16數據集來訓練我們的UNET模型。

在LUNA16數據集中，每個CT掃描用結節點和用於生成二進制掩模的結節的半徑來標注。 我將首先討論LUNA16數據集的預處理。 在數據集中，CT掃描保存在“.mhd”文件中，SimpleITK用於讀取圖像。 我已經定義了三個功能：