U-net实现(keras)

毕设老师给的题目是基于深度学习的肝脏肿瘤分割，而unet则是深度网络实现图像分割的benchmark。本人小白一名，没有机器学习与深度学习的基础，但因为毕设紧迫。来不及系统学习，博客里记录学习过程，可能有很多地方理解并不正确，但会一边学习一边更正之前的内容。

unet结构说明
其实unet就是encoder-decoder，encoder部分是重复conv-conv-maxpooling的过程，用来获取图像的高层抽象信息(此时的feature size(height & width)，就会小于原输入图像，但通道数(feature maps number)会变多，通道数越多，提取到的特征越多），decoder部分则需要还原分辨率(feature size)，是通过transpose convolution实现的，同时与上层的feature map进行concatenate（更新：这可以理解为long skip connection，unet的三维实现：vnet，它对unet进行的改进之一就是加入了short skip connection，就是我们熟悉的residual block），重复transconv-concate-conv-conv的过程，最后得到与原输入一样大小的分割结果图。

函数比较
下面两个用keras实现的代码（均来自github上的开源项目），可以发现其一使用了Conv2DTranspose,这与keras里的Deconvolution2D效果一样，也与tf.nn.conv2d_transpose效果一样。另一个则使用了UpSampling2D，虽然这两个都可以增加分辨率，但实现原理是不同的，并且，Conv2DTranspose在增加分辨率的同时，可以自定义输出通道数，在unet中，这儿的输出通道一般是减半，而UpSampling2D仅能增加分辨率，不能改变输出通道数。（第一个实现与原unet论文提出的结构更符合）

keras对比tensorflow
不知道自己的理解是否正确，但是经过比较发现，对于tensorflow建立的model,一般只到activation,模型的loss与optimizer都是另外定义的（有时候封装成一个函数）。而keras的模型，一般直接建立到底，即在activation后，会跟着model.compile(optimizer,loss,metrics,...)，猜测，应该是tensorflow需要sess.run([需要的结果],feed_dict)，对于训练和测试，两个的feed_dict不同，有时需要的结果也不同，因此模型的结构、损失与优化最好分模块写，而keras封装了model.fit用于训练，model.predict用于预测，所以直接把整个模型结构、损失与优化写在一起。
更新：其实，两者都是类似的，都需要单独的几个东西：
一是原始的y_pred值，这个值在本案例中进行了sigmoid，有些网络中没有进行激活，若没有激活的话，那么后期选择的损失函数与评价函数，就得加上激活的那一部分。
二是损失函数，即优化器优化的对象
三是评价函数，这个可以直接由损失函数取反，也可以定义不同的评价标准
keras与tensorflow只是在实现的细节上有区别。
比如keras的
model.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss, metrics=[dice_coef])
就类似于tensorflow的

loss,  acc, _ = sess.run([dice_coef_loss, dice_coef, optimizer], feed_dict={#训练数据集
}

keras的model.predict
就类似于tensorflow的

loss,  acc = sess.run([dice_coef_loss, dice_coef], feed_dict={# 测试数据集
}

学习感想
很多时候，感觉学习的难点不在于如何建立网络。而是整个的步骤，数据获取，预处理，损失函数和metrics怎么根据自己的数据定义等等。特别是学习完别人的例子后，换一个数据来源就不知道怎么处理了，也不知道如何检测自己实现的正确性，网络上完整的例子，进阶式的教学很少，数据集的获取也是个小问题，翻墙下载太慢等等。日后会继续更新博客，争取实现从零到完整的项目。

keras 实现

# from kaggle nerve segmentation competition 
def dice_coef(y_true, y_pred):
    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_unet():
    inputs = Input((img_rows, img_cols, 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)

    up6 = concatenate([Conv2DTranspose(256, (2, 2), strides=(2, 2), padding='same')(conv5), conv4], axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)

    up7 = concatenate([Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(conv6), conv3], axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)

    up8 = concatenate([Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv7), conv2], axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)

    up9 = concatenate([Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same')(conv8), 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

# from retina segmentation code
def get_unet(n_ch,patch_height,patch_width):
    inputs = Input(shape=(patch_height,patch_width,n_ch))
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(inputs)
    conv1 = Dropout(0.2)(conv1)
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv1)
    pool1 = MaxPooling2D((2, 2))(conv1)
    #
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = Dropout(0.2)(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    pool2 = MaxPooling2D((2, 2))(conv2)
    #
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = Dropout(0.2)(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)

    up1 = UpSampling2D(size=(2, 2))(conv3)
    up1 = concatenate([conv2,up1],axis=3)
    conv4 = Conv2D(64, (3, 3), activation='relu', padding='same')(up1)
    conv4 = Dropout(0.2)(conv4)
    conv4 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv4)
    #
    up2 = UpSampling2D(size=(2, 2))(conv4)
    up2 = concatenate([conv1,up2], axis=3)
    conv5 = Conv2D(32, (3, 3), activation='relu', padding='same')(up2)
    conv5 = Dropout(0.2)(conv5)
    conv5 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv5)
    #
    conv6 = Conv2D(2, (1, 1), activation='relu',padding='same')(conv5)
    # 这一步是针对该数据集的，unet原始网络中没有这个
    #conv6 = core.Reshape((2,patch_height*patch_width))(conv6)
    #conv6 = core.Permute((2,1))(conv6)
    # 经过评论区指正，conv6 = Conv2D(2, (1, 1)）输出的tensor为(n,h,w,2)
    # 原参考案例是theno后端，才会写成core.Reshape((2,patch_height*patch_width))
    conv6 = core.Reshape((-1,2))(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

后续说明
在学习的时候对于网络loss和metrics的选取有疑问，这两个的输入都是灰度图，第二个例子对于mask（即groundtruth）进行了变形，将(N,h,w,1) -->(N,hw,2)，最后一个维度变成了2是对class进行了one-hot编码，即前景编码为01，背景编码为10。这样才可以使用categorical_crossentropy，与对应的metrics=[‘accuracy’]。
查找资料，stackoverflow上说，对于multiclass的分类，有几个class，最后就需要对应几个feature map（即channel数量），一个channel对应一个class的mask，1代表为该class，0代表是其他的class，并使用loss=‘categorical_crossentropy’,metrics=[‘accuracy’]。
举例：原始mask为(100,64,64),100张，64的宽高。
其中0代表背景，1代表汽车，2代表人。那么在提供给网络时，
应先mask = to_categorical(mask, 3),
mask的shape变为(10064*64, 3)，
再mask = mask.reshape((100, 64, 64, 3))