深度學習（七）U-Net原理以及keras程式碼實現醫學影象眼球血管分割
原文作者：aircraft
原文連結：https://www.cnblogs.com/DOMLX/p/9780786.html
DRIVE資料集下載百度雲連結：連結：https://pan.baidu.com/s/1C_1ikDwexB0hZvOwMSeDtw
提取碼：8m1q
U-net+kears實現眼部血管分割原始碼python2.7版本的百度雲連結： 連結：https://pan.baidu.com/s/1C_1ikDwexB0hZvOwMSeDtw
提取碼：8m1q
U-net+kears實現眼部血管分割原始碼python3.6版本的百度 連結：連結：https://pan.baidu.com/s/1rAf6wuWGCswuBfkDivxjyQ
提取碼：rtdg
全卷積神經網路
大名鼎鼎的FCN就不多做介紹了，這裡有一篇很好的博文 ofollow,noindex" target="_blank">http://www.cnblogs.com/gujianhan/p/6030639.html 。
不過還是建議把論文讀一下，這樣才能加深理解。
醫學影象分割框架
醫學影象分割主要有兩種框架，一個是基於CNN的，另一個就是基於FCN的。這裡都是通過網路來進行語義分割。
那麼什麼是語義分割？可不是漢字分割句意，在影象處理中有自己的定義。
影象語義分割的意思就是機器自動分割並識別出影象中的內容，比如給出一個人騎摩托車的照片，機器判斷後應當能夠生成右側圖，紅色標註為人，綠色是車（黑色表示 back ground）。
所以影象分割對影象理解的意義， 就好比讀古書首先要斷句一樣 。
在 Deeplearning 技術快速發展之前，就已經有了很多做影象分割的技術，其中比較著名的是一種叫做 “Normalized cut” 的圖劃分方法，簡稱 “ N-cut ”。
N-cut 的計算有一些連線權重的公式，這裡就不提了，它的思想主要是 通過畫素和畫素之間的關係權重來綜合考慮，根據給出的閾值，將影象一分為二 。
基於CNN 的框架
這個想法也很簡單，就是對影象的每一個畫素點進行分類，在每一個畫素點上取一個patch，當做一幅影象，輸入神經網路進行訓練，舉個例子：
這是一篇發表在NIPS上的論文 Ciresan D, Giusti A, Gambardella L M, et al. Deep neural networks segment neuronal membranes in electron microscopy images[C]//Advances in neural information processing systems. 2012: 2843-2851.
這是一個二分類問題，把影象中所有label為0的點作為負樣本，所有label為1的點作為正樣本。
這種網路顯然有兩個缺點：
1. 冗餘太大，由於每個畫素點都需要取一個patch，那麼相鄰的兩個畫素點的patch相似度是非常高的，這就導致了非常多的冗餘，導致網路訓練很慢。
2. 感受野和定位精度不可兼得，當感受野選取比較大的時候，後面對應的pooling層的降維倍數就會增大，這樣就會導致定位精度降低，但是如果感受野比較小，那麼分類精度就會降低。
基於FCN框架
在醫學影象處理領域，有一個應用很廣泛的網路結構—- U-net ,網路結構如下：
可以看出來，就是一個全卷積神經網路，輸入和輸出都是影象，沒有全連線層。較淺的高解析度層用來解決畫素定位的問題，較深的層用來解決畫素分類的問題。
好了理解完U-net網路，我們就學習一下怎麼用U-net網路來進行醫學影象分割。
U-net+kears實現眼部血管分割
原作者的【 英文說明 】https://github.com/orobix/retina-unet#retina-blood-vessel-segmentation-with-a-convolution-neural-network-u-net
實現環境可直接看這篇部落格下載： 2018最新win10 安裝tensorflow1.4（GPU/CPU）+cuda8.0+cudnn8.0-v6 + keras 安裝CUDA失敗 匯入tensorflow失敗報錯問題解決
linux下就環境一樣，配置就要自己去找了。
1、介紹
為了能夠更好的對眼部血管等進行檢測、分類等操作，我們首先要做的就是對眼底影象中的血管進行分割，保證最大限度的分割出眼部的血管。從而方便後續對血管部分的操作。
這部分程式碼選用的資料集是DRIVE資料集，包括訓練集和測試集兩部分。眼底影象資料如圖1所示。
圖1 DRIVE資料集的訓練集眼底影象
DRIVE資料集的優點是：不僅有已經手工分好的的血管影象（在manual資料夾下，如圖2所示），而且還包含有眼部輪廓的影象（在mask資料夾下，如圖3所示）。
圖2 DRIVE資料集的訓練集手工標註血管影象
圖3 DRIVE資料集的訓練集眼部輪廓影象
DRIVE資料集的缺點是：顯而易見，從上面的圖片中可以看出，訓練集只有20幅圖片，可見資料量實在是少之又少。。。
所以，為了得到更好的分割效果，我們需要對這20幅影象進行預處理從而增大其資料量
2、依賴的庫
- 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
3、資料讀取與儲存
資料集中訓練集和測試集各只有20幅眼底影象（tif格式）。首先要做的第一步就是對生成資料檔案，方便後續的處理。所以這裡我們需要對資料集中的眼底影象、人工標註的血管影象、眼部輪廓生成資料檔案。這裡使用的是hdf5檔案。有關hdf5檔案的介紹，請參考CSDN部落格（ HDF5快速上手全攻略 ）。
4、網路解析
因為U-net網路可以針對很少的資料集來進行語義分割，比如我們這個眼球血管分割就是用了20張圖片來訓練就可以達到很好的效果。而且我們這種眼球血管，或者指靜脈，指紋之類的提取特徵或者血管靜脈在U-net網路裡就是一個二分類問題，大家一聽，二分類對於目前的神經網路不是一件很簡單的事情了嗎？還有是什麼可以說的。
的確目前二分類問題是沒有什麼難度了，只要給我足夠的資料集做訓練。而本文用的U-net網路來實現這個二分類就只需要二十張圖片來作為資料集。大家可以看到優勢所在了吧。
5、具體實現
首先我們肯定都是要對資料進行一些預處理的。
第一步
先將影象轉為灰度圖分別讀入陣列建立起一個符合我們自己的tensor的格式才好傳入神經網路，這裡我們是先將資料存入 hdf5檔案 中，在開始執行的時候從檔案中讀入。
#將對應的影象資料存入對應影象陣列
def get_datasets(imgs_dir,groundTruth_dir,borderMasks_dir,train_test="null"):
imgs = np.empty((Nimgs,height,width,channels))
groundTruth = np.empty((Nimgs,height,width))
border_masks = np.empty((Nimgs,height,width))
for path, subdirs, files in os.walk(imgs_dir): #list all files, directories in the path
for i in range(len(files)):
#original
print("original image: " +files[i])
img = Image.open(imgs_dir+files[i])
imgs[i] = np.asarray(img)
#corresponding ground truth
groundTruth_name = files[i][0:2] + "_manual1.gif"
print("ground truth name: " + groundTruth_name)
g_truth = Image.open(groundTruth_dir + groundTruth_name)
groundTruth[i] = np.asarray(g_truth)
#corresponding border masks
border_masks_name = ""
if train_test=="train":
border_masks_name = files[i][0:2] + "_training_mask.gif"
elif train_test=="test":
border_masks_name = files[i][0:2] + "_test_mask.gif"
else:
print("specify if train or test!!")
exit()
print("border masks name: " + border_masks_name)
b_mask = Image.open(borderMasks_dir + border_masks_name)
border_masks[i] = np.asarray(b_mask)
print("imgs max: " +str(np.max(imgs)))
print("imgs min: " +str(np.min(imgs)))
assert(np.max(groundTruth)==255 and np.max(border_masks)==255)
assert(np.min(groundTruth)==0 and np.min(border_masks)==0)
print("ground truth and border masks are correctly withih pixel value range 0-255 (black-white)")
#reshaping for my standard tensors
imgs = np.transpose(imgs,(0,3,1,2))
assert(imgs.shape == (Nimgs,channels,height,width))
groundTruth = np.reshape(groundTruth,(Nimgs,1,height,width))
border_masks = np.reshape(border_masks,(Nimgs,1,height,width))
assert(groundTruth.shape == (Nimgs,1,height,width))
assert(border_masks.shape == (Nimgs,1,height,width))
return imgs, groundTruth, border_masks
if not os.path.exists(dataset_path):
os.makedirs(dataset_path)
#getting the training datasets
imgs_train, groundTruth_train, border_masks_train = get_datasets(original_imgs_train,groundTruth_imgs_train,borderMasks_imgs_train,"train")
print("saving train datasets")
write_hdf5(imgs_train, dataset_path + "DRIVE_dataset_imgs_train.hdf5")
write_hdf5(groundTruth_train, dataset_path + "DRIVE_dataset_groundTruth_train.hdf5")
write_hdf5(border_masks_train,dataset_path + "DRIVE_dataset_borderMasks_train.hdf5")
#getting the testing datasets
imgs_test, groundTruth_test, border_masks_test = get_datasets(original_imgs_test,groundTruth_imgs_test,borderMasks_imgs_test,"test")
print("saving test datasets")
write_hdf5(imgs_test,dataset_path + "DRIVE_dataset_imgs_test.hdf5")
write_hdf5(groundTruth_test, dataset_path + "DRIVE_dataset_groundTruth_test.hdf5")
write_hdf5(border_masks_test,dataset_path + "DRIVE_dataset_borderMasks_test.hdf5")
第二步
是對讀入記憶體準備開始訓練的影象資料進行一些增強之類的處理，這裡對其進行了，直方圖均衡化，資料標準化，並且壓縮畫素值到0-1，將其的一個數據符合標準正態分佈。當然啦我們這個資料拿來訓練還是太少的，所以我們對每張圖片取patch時，除了正常的每個 patch 每個 patch移動的取之外，我們還在資料範圍內進行隨機取 patch ，這樣雖然各個 patch 之間會有一部分資料是相同的，但是這對於網路而言，你傳入的也是一個新的東西，網路能從中提取到的特徵也更多了。這一步的目的其實就是在有限的資料集中進行一些資料擴充，這也是在神經網路訓練中常用的手段了。
當然了在這個過程中我們也可以隨機組合小的patch來看看。
隨機原圖：
mask圖：
處理待訓練資料的部分程式碼：
def get_data_training(DRIVE_train_imgs_original,
DRIVE_train_groudTruth,
patch_height,
patch_width,
N_subimgs,
inside_FOV):
train_imgs_original = load_hdf5(DRIVE_train_imgs_original)
train_masks = load_hdf5(DRIVE_train_groudTruth) #masks always the same
# visualize(group_images(train_imgs_original[0:20,:,:,:],5),'imgs_train')#.show()#check original imgs train
train_imgs = my_PreProc(train_imgs_original)#直方圖均衡化，資料標準化，壓縮畫素值到0-1
train_masks = train_masks/255.
train_imgs = train_imgs[:,:,9:574,:]#cut bottom and top so now it is 565*565
train_masks = train_masks[:,:,9:574,:]#cut bottom and top so now it is 565*565
data_consistency_check(train_imgs,train_masks)
#check masks are within 0-1
assert(np.min(train_masks)==0 and np.max(train_masks)==1)
print("\ntrain images/masks shape:")
print(train_imgs.shape)
print("train images range (min-max): " +str(np.min(train_imgs)) +' - '+str(np.max(train_imgs)))
print("train masks are within 0-1\n")
#extract the TRAINING patches from the full images
patches_imgs_train, patches_masks_train = extract_random(train_imgs,train_masks,patch_height,patch_width,N_subimgs,inside_FOV)
data_consistency_check(patches_imgs_train, patches_masks_train)
print("\ntrain PATCHES images/masks shape:")
print(patches_imgs_train.shape)
print("train PATCHES images range (min-max): " +str(np.min(patches_imgs_train)) +' - '+str(np.max(patches_imgs_train)))
return patches_imgs_train, patches_masks_train#, patches_imgs_test, patches_masks_test
第三步
按照U-net的網路結構，使用keras來構造出網路。這裡對keras函式語法不太理解的可以看這篇部落格： 深度學習（六）keras常用函式學習
def get_unet(n_ch,patch_height,patch_width):
inputs = Input(shape=(n_ch,patch_height,patch_width))
conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(inputs)
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)
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)
conv6 = core.Permute((2,1))(conv6)
############
conv7 = core.Activation('softmax')(conv6)
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
第四步
這裡就是將我們處理好的資料傳入到網路裡訓練了。得出結果圖。
看的出來很多很細的紋理都被提取出來了，這個U-net網路也可以用於一些醫學細胞的 邊緣提取，指靜脈，掌靜脈之類的紋路提取都可以。後面在下可能還會出指靜脈之類其他影象的語義分割提取，關注在下就可以看到啦hhhhhhh
因為keras內部可以直接將整個網路結構打印出來，所以我們可以看到完整的網路結構圖如下所示：
參考部落格：https://blog.csdn.net/u013063099/article/details/79981097
參考部落格：https://blog.csdn.net/qq_16900751/article/details/78251778