論文復現之醫學影像應用：視網膜血管分割

2019 年 10 月 6 日
筆記

論文復現之醫學影像應用：視網膜血管分割

0.導語

今日研究為繼續上次論文中的一個內容：U-Net網路，於是找了一篇經典論文，並學習論文及程式碼解讀。在學習U-Net網路後，使用U-Net神經網路提取視網膜紋理血管。

1.論文閱讀

論文題目：U-Net Convolutional Networks for Biomedical Image Segmentation。

中文翻譯為：用於生物醫學影像分割的U-Net卷積網路。

1.1 摘要

之前，在訓練一個深度網路需要大量已標註的訓練樣本。在這篇論文中，提出了一種網路和訓練策略。為了更有效的使用標註數據，使用了數據擴張的方法。本片論文的網路分為兩部分：收縮路徑(contracting path)和擴張路徑(expanding path)。

收縮路徑主要用於捕捉上下文特徵
擴展路徑用於定位

該網路獲得了贏得了ISBI cell tracking challenge 2015。不僅如此，這個網路非常的快，對一個512×512的影像，使用一塊GPU只需要不到一秒的時間。

1.2 介紹

卷積神經網路的典型用途是分類任務，其中輸出到影像是單個類別標籤。然而，在許多視覺任務中，尤其是在生物醫學影像處理中，期望的輸出應該包括定位，即：應該將類別標籤分配給每個像素。此外，在生物醫學任務中，千量級的訓練影像通常難以訓練。因此，Ciresan等人在滑動窗口設置中訓練網路，通過提供圍繞該像素的局部區域（修補程式）作為輸入來預測每個像素的類別標籤。

首先，這個網路可以局部化。
其次，修補程式方面的訓練數據遠大於訓練影像的數量。

由此產生的網路大幅度贏得了ISBI 2012的EM segmentation challenge。

顯然，Ciresan等人的策略有兩個缺點。

首先，它非常慢。

因為網路必須分別為每個修補程式運行，並且由於修補程式重疊而導致大量冗餘。

其次，需要在局部準確性和獲取整體上下文資訊之間取捨。

較大的修補程式需要更多的最大池化層來降低局部準確性，而較小的修補程式則使網路只能看到很少的上下文。最近的方法提出了一個分類器輸出，它考慮了來自多個層的特徵。 良好的本地化和上下文的使用是可能的同時。

在這篇文章中，就建立了一個更加優雅的框架，通常被稱為「全卷積網路」（FCN)。隨後修改並拓展了這個框架，使其可以僅使用少量訓練圖片就可以工作，獲得更高的分割準確率。網路架構如下圖所示：

對FCN中的重要修改是：

上取樣部分有大量的特徵通道，它們允許網路將上下文資訊傳播到更高解析度的層。

使得網路架構中擴張路徑與收縮路徑基本對稱，併產生u形結構。

沒有任何完全連接的層，分割圖僅包含像素，對於該像素，輸入影像中的完整上下文是可用的。

該策略允許通過重疊拼貼策略對任意大的影像進行無縫分割。 為了預測影像的邊界區域中的像素，通過鏡像輸入影像來推斷丟失的上下文。這種平鋪策略對於將網路應用於大影像很重要，否則解析度將受到GPU記憶體的限制。

在上述談了，使用少量數據做訓練，在這篇論文中也使用了數據增強方法，對於這篇論文中的任務，論文中通過對可用的訓練影像應用彈性變形來使用過多的數據增強。

這使得網路能夠學習這種變形的不變性，而不需要在注釋影像語料庫中看到這些變換。這在生物醫學分割中特別重要，因為變形曾經是組織中最常見的變異，並且可以有效地模擬真實變形。 Dosovitskiy等人證明了學習不變性的數據增強在無監督特徵學習的範圍內的價值。

許多細胞分割任務中的另一個挑戰是分離相同類別的觸摸物體。

為此，論文中提出使用加權損失：觸摸單元之間的分離背景標籤在損失函數中獲得大的權重。

1.3 網路架構

網路架構如上圖所示：它由一條收縮路徑（左側）和一條擴展路徑（右側）組成。

【收縮路徑】

收縮路徑是典型的卷積網路架構。它的架構是一種重複結構，每次重複中都有2個卷積層和一個pooling層，卷積層中卷積核大小均為3×3，激活函數使用ReLU，兩個卷積層之後是一個2×2的步長為2的最大池化層。每一次下取樣後我們都把特徵通道的數量加倍。

【擴張路徑】

擴展路徑中的每一步包括對特徵映射先進行上取樣，然後進行2×2卷積（「上卷積」）。該特徵映射將特徵通道的數量減半，與收縮路徑中相應裁剪的特徵拼接起來，以及兩個3×3卷積，每個卷積都有一個ReLU。由於每個卷積中邊界像素的丟失，裁剪是必要的。在最後一層，使用1×1卷積來將每個64通道特徵圖映射到特定深度的結果。

【卷積層總數】

網路總共有23個卷積層。

1.4 訓練

論文中採用隨機梯度下降法(SGD)訓練。由於沒有使用0填補的卷積，輸出影像比輸入小一個恆定的邊界寬度。為了最大限度地降低開銷並最大限度地利用GPU記憶體，我們傾向於在較大批量的情況下使用較大的輸入切片，從而將批量減少為單個影像。因此，我們使用高動量（0.99），以便大量先前看到的訓練樣本確定當前優化步驟中的更新。

損失函數使用逐像素的 softmax 函數和交叉熵損失函數的結合
Softmax函數：

a_k(x)表示在feature maps中的的channel=k的feature map像素位置為x的激活值。K表示總類別數，pk(x)表示似然函數。

交叉熵函數：

其中l是每個像素的真實標籤，w是權重地圖，表示訓練中某些像素更加重要。

為了使某些像素點更加重要，我們在公式中引入了w(x)。我們對每一張標註影像預計算了一個權重圖，來補償訓練集中每類像素的不同頻率，使網路更注重學習相互接觸的細胞之間的小的分割邊界。我們使用形態學操作計算分割邊界。權重圖計算公式如下：

其中w_c是權重地圖來平衡類像素的頻率。d1表示最近單元邊界的距離。d2表示到第二最近單元的邊界的距離。文中設置w_0=10,σ≈5pixels。

訓練需要標註對應的mask,就是類別的區域標記。

網路中權重的初始化：在具有許多卷積層和通過網路的不同路徑的深層網路中，權重的良好初始化非常重要。否則，網路的某些部分可能會過度激活，而其他部分則無法提供。理想情況下，應調整初始權值，使網路中的每個特徵映射具有近似單位差異。對於具有我們的體系結構（交替卷積和ReLU層）的網路，這可以通過從具有標準偏差(2/N)^0.5的高斯分布中繪製初始權重來實現，其中N為每個神經元的輸入節點數量。例如3×3的64通道的卷積層N=3x3x64=576。

1.5 數據增加

在只有少量樣本的情況況下，要想儘可能的讓網路獲得不變性和魯棒性，數據增強是必不可少的。在顯微影像的情況下，我們主要需要移位和旋轉不變性以及對變形和灰度值變化的魯棒性。論文中使用隨機位移矢量在粗糙的3×3網格上產生平滑形變。位移是從10像素標準偏差的高斯分布中取樣的。然後使用雙三次插值計算每個像素的位移。在壓縮路徑末端的退出層執行進一步的隱式數據增強。

最後論文的實驗部分，這裡直接在DRIVE資料庫上做實驗！

2.視網膜血管分割實驗

2.1 實驗任務

實驗任務：使用U-Net神經網路提取紋理血管。

為什麼要做這個，有什麼實際意義？

臨床實驗中我們要能夠更好的對眼部血管等進行檢測、分類等操作，我們首先要做的就是對眼底影像中的血管進行分割，保證最大限度的分割出眼部的血管。從而方便後續對血管部分的操作。

2.2 數據集

實驗數據集為：DRIVE數據集。

數據集下載地址： http://www.isi.uu.nl/Research/Databases/DRIVE/

【數據集識別】

數據集分為訓練集與測試集，而每個都分為如下目錄：

1st_manual目錄下的數據集為：手工分好的血管影像，如下圖所示：

images目錄下的數據集為眼底影像，如下圖所示：

mask目錄下的數據集為眼部輪廓影像，如下圖所示：

從上面可以看出數據集非常少，總共只有20幅圖片，為此我們需要做預處理增大數據集。

【U-Net網路優勢】

在上述U-net論文中提到U-Net網路可以針對很少的數據集來進行語義分割，比如我們這個眼球血管分割就是用了20張圖片來訓練就可以達到很好的效果。而且我們這種眼球血管，或者指靜脈，指紋之類的提取特徵或者血管靜脈在U-net網路里就是一個二分類問題。而本文用的U-net網路來實現這個二分類就只需要二十張圖片來作為數據集。

2.2 運行環境及庫

復現論文程式碼：

官方：https://github.com/orobix/retina-unet 非官方：https://github.com/CNU105/Retina-Unet

numpy >= 1.11.1
keras >= 2.1.0
PIL >=1.1.7
opencv >=2.4.10
h5py >=2.6.0
configparser >=3.5.0b2
scikit-learn >= 0.17.1

上面是官方給出的運行庫環境。實際上當你安裝上後還是有問題的！

還需要安裝如下內容：

graphviz
pydot
pydot-ng

如果graphviz安裝報錯，則需要：

sudo apt-get install graphviz或者yum安裝

但是我的是伺服器上，沒有root許可權，怎麼安裝呢？

這就需要源碼編譯，在linux系統中，源碼編譯基本可以解決一切安裝問題！

由於沒有安裝graphviz，同時需要root許可權，所以安裝就出現了問題，然後只能通過源碼編譯，結果源碼編譯最新版也出現了問題，說是最新版的bug，我就轉向系統的yum本地用戶安裝：

# 查看yum中是否有你需要的包  yum list 'graphviz*'  # 下載rpm包(上述輸出：graphviz.x86_64)  yumdownloader graphviz.x86_64  # 解壓RPM包  rpm2cpio graphviz-2.30.1-19.el7.x86_64.rpm |cpio -idvm  # 添加環境變數  vim  ~/.bashrc  export PATH=$PATH:$HOME/usr/bin/  source ~/.bashrc  # 備註：這種安裝會在當前目錄下新建一個usr目錄，安裝的軟體全在usr下的bin目錄下！

安裝好上述就可以了。

但是我的伺服器上沒python，我也沒root許可權，如何安裝python呢?

那就是上述說的源碼編譯，指定源碼編譯路徑，然後添加環境變數即可！

過程中重要兩點：

（1）編譯路徑指定

./configure --prefix=/home/city/python3/

編譯後，直接make，然後make install即可。

（2）配置環境變數

vim  ~/.bashrc  export PATH=$PATH:$HOME/home/city/python3/  source ~/.bashrc

至此，通過源碼編譯與本地用戶安裝解決了環境問題！

2.3 開始實驗

【數據預處理】

上面說數據集的時候，提到數據集非常少，那麼在實驗前就得做數據預處理，而該數據預處理做的工作就是：讀取原來的訓練集與測試集，重新生成數據文件，並保存為hdf5文件。

對應到程式碼中就是運行下面程式碼：

python3 prepare_datasets_DRIVE.py

運行該程式碼後，會生成如下目錄及文件：

【配置文件】

在程式碼中有個configuration.txt配置文件，該配置文件主要用來配置數據集目錄，訓練及測試設置！

【開始訓練】

運行程式碼

在官網中是讓我們運行下面程式碼：

python3 run_training.py

而實際上直接可以運行：

python3 ./src/retinaNN_training.py

因為，上述程式碼實際上就是獲取我們的實驗設置參數，判斷是否是GPU訓練而已。如果自己會後台運行，直接敲命令即可！

通過運行retinaNN_training.py我們的實驗便開始進入訓練過程了，訓練過程非常漫長，本次運行GPU共8塊，運行時長達到6個h。

模型查看

網路定義每行這個小括弧填inputs是代表這層模型連接在inputs之後

inputs = Input(shape=(n_ch,patch_height,patch_width))  # 2D 卷積層  conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(inputs)

完整模型注釋及程式碼：

def get_unet(n_ch,patch_height,patch_width):      # 第一層      inputs = Input(shape=(n_ch,patch_height,patch_width))      # 2D 卷積層      conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(inputs)      # Dropout拋棄一些參數防止過擬合      conv1 = Dropout(0.2)(conv1)      conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv1)      pool1 = MaxPooling2D((2, 2))(conv1)        # 第二層      conv2 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(pool1)      conv2 = Dropout(0.2)(conv2)      conv2 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv2)      pool2 = MaxPooling2D((2, 2))(conv2)        # 第三層      conv3 = Conv2D(128, (3, 3), activation='relu', padding='same',data_format='channels_first')(pool2)      conv3 = Dropout(0.2)(conv3)      conv3 = Conv2D(128, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv3)      # 上取樣函數，size(2,2)其實就等於將原圖放大四倍（水平兩倍，垂直兩倍）例如32*32 變成 64*64的影像      up1 = UpSampling2D(size=(2, 2))(conv3)        # 兩層拼接      up1 = concatenate([conv2,up1],axis=1)        conv4 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(up1)      conv4 = Dropout(0.2)(conv4)      conv4 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv4)      up2 = UpSampling2D(size=(2, 2))(conv4)        # # 兩層拼接      up2 = concatenate([conv1,up2], axis=1)        conv5 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(up2)      conv5 = Dropout(0.2)(conv5)      conv5 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv5)      conv6 = Conv2D(2, (1, 1), activation='relu',padding='same',data_format='channels_first')(conv5)      conv6 = core.Reshape((2,patch_height*patch_width))(conv6)      # Permute層將輸入的維度按照給定模式進行重排      conv6 = core.Permute((2,1))(conv6)      conv7 = core.Activation('softmax')(conv6)        # 將模型的輸入和輸出給Model函數就會自己組建模型運行圖結構      model = Model(inputs=inputs, outputs=conv7)      # sgd = SGD(lr=0.01, decay=1e-6, momentum=0.3, nesterov=False)      # 配置模型      model.compile(optimizer='sgd', loss='categorical_crossentropy',metrics=['accuracy'])      return model

模型訓練程式碼流程：

首先定義模型
讀取配置文件
載入數據

這一步非常重要，對讀入記憶體準備開始訓練的影像數據進行一些增強之類的處理。這裡對其進行了，直方圖均衡化，數據標準化，並且壓縮像素值到0-1，將其的一個數據符合標準正態分布。(這個處理程式碼在lib目錄下的extract_patches.py中)。

前面我們知道數據集非常少，那怎麼讓訓練的數據集增大呢？

這裡就是做這個處理，即對每張圖片取patch時，除了正常的每個patch每個patch移動的取之外，我們還在數據範圍內進行隨機取patch，這樣雖然各個patch之間會有一部分數據是相同的，但是這對於網路而言，你傳入的也是一個新的東西，網路能從中提取到的特徵也更多了。這一步的目的其實就是在有限的數據集中進行一些數據擴充，這也是在神經網路訓練中常用的手段了。

存儲隨機組合的patch

左圖為隨機原圖，右圖為mask圖

存儲模型

【測試及評估】

python3 src/retinaNN_predict.py

由於伺服器上沒有可視化頁面，所以當調用matlibplot時候，會報如下錯誤，那麼就需要後台運行，不顯示！

no display name and no $DISPLAY environment variable

解決辦法是：

加入下面兩行程式碼即可！

import matplotlib  matplotlib.use('Agg')

在本程式碼中，做了如下工作：

同訓練集一樣取patch
利用訓練好的模型載入測試數據
存儲並可視化結果
評估結果

從左到右，依次為真實，原圖，預測。

從上到下，依次為真實，原圖，預測。