在之前的一篇博客里：医学图像分割 unet实现(一)，是学习并复现别人的实验。这篇将记录下自己毕设第一阶段的实验。毕设题目为：基于深度学习的肝脏肿瘤分割。
经过几番调整，最终确定：第一阶段分割出腹部图像中的肝脏，作为第二阶段的ROI(region of interest)，第二阶段利用ROI对腹部图像进行裁剪，裁剪后的非ROI区域变成黑色，作为该阶段输入，分割出肝脏中的肿瘤（更新2019-4-2，已做完实验，医学图像分割 基于深度学习的肝脏肿瘤分割 实战(二）。第三阶段用随机场的后处理方法进行优化。
此外，觉得当时的那篇文章写得很不用心，没费多大功夫，复制粘贴为主。我记录的主要原因之一就是觉得这方面入门实战的文章很少，希望让情况类似的同学少走一点弯路。可现在自己写的东西过几天再去看，都会一头雾水。没有从主干思路到细节的分层讲解，没有关键点的介绍。这篇文章开始，将尽我所能，写得利于接受一点。
（不过，自己才开始学习一个多月，python和框架都是临时学的，并不能保证博客里的东西是对的，有错误的话，望指出）

正文：

目标

分割出CT腹部图像的肝脏区域。

原始数据介绍

实验用到的数据为3Dircadb，即腹部的CT图。一个病人为一个文件夹。
例如3Dircadb1.1（一共20位），该文件夹下会使用到的子文件夹为PATIENT_DICOM（病人原始腹部3D CT图），MASKS_DICOM（该文件夹下具有不同部位的分割结果，比如liver和liver tumor等等）。如下图所示：

PATIENT_DICOM利用软件展示效果如下：一个dcm文件包含129张切片。

MASKS_DICOM下的liver分割图效果如下：

整体思路

1、数据提取

数据读取：
从原始dcm格式读入成我们需要的数组格式

#part1
import numpy as np
import pydicom
import os
import matplotlib.pyplot as plt
import cv2
from keras.preprocessing.image import ImageDataGenerator
from HDF5DatasetWriter import HDF5DatasetWriter
from HDF5DatasetGenerator import HDF5DatasetGenerator

for i in range(1,18): # 前17个人作为测试集
   full_images = [] # 后面用来存储目标切片的列表
   full_livers = [] #功能同上
   # 注意不同的系统，文件分割符的区别
   label_path = '～/3Dircadb/3Dircadb1.%d/MASKS_DICOM/liver'%i
   data_path = '～/3Dircadb/3Dircadb1.%d/PATIENT_DICOM'%i
   liver_slices = [pydicom.dcmread(label_path + '/' + s) for s in os.listdir(label_path)]
   # 注意需要排序，即使文件夹中显示的是有序的，读进来后就是随机的了
   liver_slices.sort(key = lambda x: int(x.InstanceNumber))
   # s.pixel_array 获取dicom格式中的像素值
   livers = np.stack([s.pixel_array for s in liver_slices])
   image_slices = [pydicom.dcmread(data_path + '/' + s) for s in os.listdir(data_path)]
   image_slices.sort(key = lambda x: int(x.InstanceNumber))
   
   """ 省略进行的预处理操作，具体见part2"""
   
   full_images.append(images)
   full_livers.append(livers)
   
   full_images = np.vstack(full_images)
   full_images = np.expand_dims(full_images,axis=-1)
   full_livers = np.vstack(full_livers)
   full_livers = np.expand_dims(full_livers,axis=-1)

2、数据的预处理

1、将ct值转化为标准的hu值
至于为什么要将值进行转化，这儿就不详细说明，具体可以参考医学图像预处理（三）——windowing(ct对比增强），这篇博文中有一样的get_pixels_hu函数
2、窗口化操作
医学图像预处理（三）——windowing(ct对比增强）
3、直方图均衡化

def clahe_equalized(imgs,start,end):
   assert (len(imgs.shape)==3)  #3D arrays
   #create a CLAHE object (Arguments are optional).
   clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))
   imgs_equalized = np.empty(imgs.shape)
   for i in range(start, end+1):
       imgs_equalized[i,:,:] = clahe.apply(np.array(imgs[i,:,:], dtype = np.uint8))
   return imgs_equalized

4、归一化
5、仅提取腹部所有切片中包含了肝脏的那些切片，其余的不要
医学图像预处理（四）—— 提取包含目标的切片

#part2
    # 接part1
   images = get_pixels_hu(image_slices)
   
   images = transform_ctdata(images,500,150)
   
   start,end = getRangImageDepth(livers)
   images = clahe_equalized(images,start,end)
   
   images /= 255.
   # 仅提取腹部所有切片中包含了肝脏的那些切片，其余的不要
  
   total = (end - 4) - (start+4) +1
   print("%d person, total slices %d"%(i,total))
   # 首和尾目标区域都太小，舍弃
   images = images[start+5:end-5]
   print("%d person, images.shape:(%d,)"%(i,images.shape[0]))
   
   livers[livers>0] = 1
   
   livers = livers[start+5:end-5]

3、数据增强

利用keras的数据增强接口，可以实现分割问题的数据增强。一般的增强是分类问题，这种情况，只需要对image变形，label保持不变。但分割问题，就需要image和mask进行同样的变形处理。具体怎么实现，参考下面代码，注意种子设定成一样的。

# 可以在part1之前设定好（即循环外）
seed=1
data_gen_args = dict(rotation_range=3,
                    width_shift_range=0.01,
                    height_shift_range=0.01,
                    shear_range=0.01,
                    zoom_range=0.01,
                    fill_mode='nearest')

image_datagen = ImageDataGenerator(**data_gen_args)
mask_datagen = ImageDataGenerator(**data_gen_args)

#part3 接part2
    image_datagen.fit(full_images, augment=True, seed=seed)
    mask_datagen.fit(full_livers, augment=True, seed=seed)
    image_generator = image_datagen.flow(full_images,seed=seed)
    mask_generator = mask_datagen.flow(full_livers,seed=seed)

    train_generator = zip(image_generator, mask_generator)
    x=[]
    y=[]
    i = 0
    for x_batch, y_batch in train_generator:
        i += 1
        x.append(x_batch)
        y.append(y_batch)
        if i>=2: # 因为我不需要太多的数据
            break
    x = np.vstack(x)
    y = np.vstack(y)

4、数据存储

一般而言，数据量较大的话，都会先将原始数据库的东西转化为np或者h5格式的文件，我感觉这样有两个好处，一是真正输入网络训练的时候io量会大大减少(特别是h5很适用于大的数据库）,二是数据分享或者上传至服务器时也方便一点。

实验中会出现两个类，分别是写h5和读h5文件的辅助类：
这读文件的类写成了generator，这样可以结合训练网络时，keras的fit_generator来使用，降低内存开销。

class HDF5DatasetWriter:
	"""用来写数据到h5文件"""
class HDF5DatasetGenerator:
    """用来读h5文件的数据"""

它们具体的实现在python h5文件的读写，因为篇幅问题，所以这儿不详述

h5文件操作需要的包：import h5py

# 可以在part1之前设定好（即循环外）
# 这儿的数量需要提前写好，感觉很不方便，但不知道怎么改，我是先跑了之前的程序，计算了一共有多少
# 张图片后再写的，但这样明显不是好的解决方案
dataset = HDF5DatasetWriter(image_dims=(2782, 512, 512, 1),
                            mask_dims=(2782, 512, 512, 1),
                            outputPath="data_train/train_liver.h5")

#part4 接part3
    dataset.add(full_images, full_livers)
    dataset.add(x, y)
    # end of lop
dataset.close()

测试数据存储的全部过程

测试数据与上面训练数据的处理过程几乎一样，但测试数据不要进行数据增强

full_images2 = []
full_livers2 = []
for i in range(18,21):#后3个人作为测试样本
    label_path = '～/3Dircadb/3Dircadb1.%d/MASKS_DICOM/liver'%i
    data_path = '～/3Dircadb/3Dircadb1.%d/PATIENT_DICOM'%i
    liver_slices = [pydicom.dcmread(label_path + '/' + s) for s in os.listdir(label_path)]
    liver_slices.sort(key = lambda x: int(x.InstanceNumber))
    livers = np.stack([s.pixel_array for s in liver_slices])
    start,end = getRangImageDepth(livers)
    total = (end - 4) - (start+4) +1
    print("%d person, total slices %d"%(i,total))
    
    image_slices = [pydicom.dcmread(data_path + '/' + s) for s in os.listdir(data_path)]
    image_slices.sort(key = lambda x: int(x.InstanceNumber))
    
    images = get_pixels_hu(image_slices)
    images = transform_ctdata(images,500,150)
    images = clahe_equalized(images,start,end)
    images /= 255.
    images = images[start+5:end-5]
    print("%d person, images.shape:(%d,)"%(i,images.shape[0]))
    livers[livers>0] = 1
    livers = livers[start+5:end-5]

    full_images2.append(images)
    full_livers2.append(livers)
    
full_images2 = np.vstack(full_images2)
full_images2 = np.expand_dims(full_images2,axis=-1)
full_livers2 = np.vstack(full_livers2)
full_livers2 = np.expand_dims(full_livers2,axis=-1)

dataset = HDF5DatasetWriter(image_dims=(full_images2.shape[0], full_images2.shape[1], full_images2.shape[2], 1),
                            mask_dims=(full_images2.shape[0], full_images2.shape[1], full_images2.shape[2], 1),
                            outputPath="data_train/val_liver.h5")


dataset.add(full_images2, full_livers2)

print("total images in val ",dataset.close())

5、构建网络

这一部分不多说，是直接用的别人写好的Unet。
唯一进行的更改，就是Crop的那部分。因为，如果输入图片的高或者宽不能等于2^n(n>=m)，m为网络收缩路径的层数。那么出现的问题就是，下采样的时候进行了取整操作，但是，后续在扩张路径又会有上采样，且上采样的结果会和收缩路径的特征图进行拼接，即long skip connection。若没有达到上面提到的要求，会导致拼接的两个特征图size不同，报错。

# partA
import os
import sys
import numpy as np
import random
import math
import tensorflow as tf
from HDF5DatasetGenerator import HDF5DatasetGenerator
from keras.models import Model
from keras.layers import Input, concatenate, Conv2D, MaxPooling2D, Conv2DTranspose,Cropping2D
from keras.optimizers import Adam
from keras.callbacks import ModelCheckpoint
from keras import backend as K
from skimage import io


K.set_image_data_format('channels_last')
    
def dice_coef(y_true, y_pred):
    smooth = 1.
    y_true_f = K.flatten(y_true)
    y_pred_f = K.flatten(y_pred)
    intersection = K.sum(y_true_f * y_pred_f)
    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)


def dice_coef_loss(y_true, y_pred):
    return -dice_coef(y_true, y_pred)

def get_crop_shape(target, refer):
        # width, the 3rd dimension
        print(target.shape)
        print(refer._keras_shape)
        cw = (target._keras_shape[2] - refer._keras_shape[2])
        assert (cw >= 0)
        if cw % 2 != 0:
            cw1, cw2 = int(cw/2), int(cw/2) + 1
        else:
            cw1, cw2 = int(cw/2), int(cw/2)
        # height, the 2nd dimension
        ch = (target._keras_shape[1] - refer._keras_shape[1])
        assert (ch >= 0)
        if ch % 2 != 0:
            ch1, ch2 = int(ch/2), int(ch/2) + 1
        else:
            ch1, ch2 = int(ch/2), int(ch/2)

        return (ch1, ch2), (cw1, cw2)

def get_unet():
    inputs = Input((IMG_HEIGHT, IMG_WIDTH , 1))
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(inputs)
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)

    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)

    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)

    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2))(conv4)

    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(pool4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)

    up_conv5 = Conv2DTranspose(256, (2, 2), strides=(2, 2), padding='same')(conv5)
    
    ch, cw = get_crop_shape(conv4, up_conv5)
    
    crop_conv4 = Cropping2D(cropping=(ch,cw), data_format="channels_last")(conv4)
    up6 = concatenate([up_conv5, crop_conv4], axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    
    up_conv6 = Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(conv6)

    ch, cw = get_crop_shape(conv3, up_conv6)
    crop_conv3 = Cropping2D(cropping=(ch,cw), data_format="channels_last")(conv3)
    
    up7 = concatenate([up_conv6, crop_conv3], axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    
    up_conv7 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv7)
    ch, cw = get_crop_shape(conv2, up_conv7)
    crop_conv2 = Cropping2D(cropping=(ch,cw), data_format="channels_last")(conv2)

    up8 = concatenate([up_conv7, crop_conv2], axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    
    up_conv8 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv8)
    ch, cw = get_crop_shape(conv1, up_conv8)
    crop_conv1 = Cropping2D(cropping=(ch,cw), data_format="channels_last")(conv1)


    up9 = concatenate([up_conv8, crop_conv1], axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)

    conv10 = Conv2D(1, (1, 1), activation='sigmoid')(conv9)

    model = Model(inputs=[inputs], outputs=[conv10])

    model.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss, metrics=[dice_coef])

    return model

6、进行训练并测试

这其中主要包括训练文件与测试文件的读取，Checkpoint回掉函数的设定，fit_generator的使用，模型预测，并对预测结果进行保存。

# partB 接partA
IMG_WIDTH = 512
IMG_HEIGHT = 512
IMG_CHANNELS = 1
TOTAL = 2782 # 总共的训练数据
TOTAL_VAL = 152 # 总共的validation数据
# part1部分储存的数据文件
outputPath = './data_train/train_liver.h5' # 训练文件
val_outputPath = './data_train/val_liver.h5'
#checkpoint_path = 'model.ckpt'
BATCH_SIZE = 8 # 根据服务器的GPU显存进行调整

class UnetModel:
    def train_and_predict(self):
        
        reader = HDF5DatasetGenerator(dbPath=outputPath,batchSize=BATCH_SIZE)
        train_iter = reader.generator()
        
        test_reader = HDF5DatasetGenerator(dbPath=val_outputPath,batchSize=BATCH_SIZE)
        test_iter = test_reader.generator()
        fixed_test_images, fixed_test_masks = test_iter.__next__()
#   
        
        model = get_unet()
        model_checkpoint = ModelCheckpoint('weights.h5', monitor='val_loss', save_best_only=True)
        # 注：感觉validation的方式写的不对，应该不是这样弄的
        model.fit_generator(train_iter,steps_per_epoch=int(TOTAL/BATCH_SIZE),verbose=1,epochs=500,shuffle=True,
                            validation_data=(fixed_test_images, fixed_test_masks),callbacks=[model_checkpoint])
#        
        reader.close()
        test_reader.close()
        
        
        print('-'*30)
        print('Loading and preprocessing test data...')
        print('-'*30)
        
        print('-'*30)
        print('Loading saved weights...')
        print('-'*30)
        model.load_weights('weights.h5')
    
        print('-'*30)
        print('Predicting masks on test data...')
        print('-'*30)
        
        
        imgs_mask_test = model.predict(fixed_test_images, verbose=1)
        np.save('imgs_mask_test.npy', imgs_mask_test)
    
        print('-' * 30)
        print('Saving predicted masks to files...')
        print('-' * 30)
        pred_dir = 'preds'
        if not os.path.exists(pred_dir):
            os.mkdir(pred_dir)
        i = 0
        
        
        for image in imgs_mask_test:
            image = (image[:, :, 0] * 255.).astype(np.uint8)
            gt = (fixed_test_masks[i,:,:,0] * 255.).astype(np.uint8)
            ini = (fixed_test_images[i,:,:,0] *255.).astype(np.uint8)
            io.imsave(os.path.join(pred_dir, str(i) + '_ini.png'), ini)
            io.imsave(os.path.join(pred_dir, str(i) + '_pred.png'), image)
            io.imsave(os.path.join(pred_dir, str(i) + '_gt.png'), gt)
            i += 1

unet = UnetModel()
unet.train_and_predict()

模型跑的过程如图。

预测结果可视化展示