为您推荐:
精华内容
最热下载
问答
  • 5.54MB weixin_39840914 2019-08-12 08:46:03
  • 5星
    1KB weixin_42666807 2021-10-03 03:15:24
  • 最近在做医学眼底血管抽取的项目,需要用到Unet,关于Unet的介绍,全网各类博客已经讲过很多,我就不再赘述了,主要是讲讲我自己训练实际的情况,和经常遇到的错误 数据集 拿到的图片一般都是以下这样的 我拿到时,...

    简介

    最近在做医学眼底血管抽取的项目,需要用到Unet,关于Unet的介绍,全网各类博客已经讲过很多,我就不再赘述了,主要是讲讲我自己训练实际的情况,和经常遇到的错误

    数据集

    拿到的图片一般都是以下这样的

    我拿到时,总共只有60张,可以说是少的可怜了,我试过直接用60张图片做训练,完全不行,train_loss都达到0.3,预测的图片黑的一片,可以肯定需要做数据增强了,由于对原图进行数据增强的时候,需要同时修改label的图片,我们可以把label的第1通道作为原图的第3通道,得到新的merge图

    这样在做数据增强的时候,就可以同时修改原图和label图了,以下是我做数据增强的关键部分代码,需要注意的是,这里的这里我对原图进行了归一化处理,mean和std不是每一次的batch产生的,而是用所有的train图像求出来的,可是原图只有60张不够,所以我另外写了一个脚本,通过数据增强获得了大概1700张图片,再求mean和std,保存成npy,以后训练和预测,都需要读取进来使用

    img = load_img(os.path.join(self.image_path, filename), \
                         target_size=(self.target_size, self.target_size))
    img = img_to_array(img)
    x_t = img.reshape((1,) + img.shape)
    img = self.datagen.flow(x_t, batch_size=1)
    batch = img.next()[0]
    train = batch[:, :, 0]
    mask = batch[:, :, 2]
    
    X[i] = img
    y[i] = mask
    
    X = X.astype('float32')
    y = y.astype('float32')
    X /= 255
    
    #X -= np.mean(x, axis=0)
    #X /= np.std(x, axis=0)
    X -= mean
    X /= std
    
    y /= 255
    y[y >= 0.5] = 1
    y[y < 0.5] = 0
    

    模型训练

    Unet模型,我是通过原模型修改过来的,原模型训练loss是0.5,acc是97,修改的模型loss是0.38,acc是98,在实际预测中抽取血管效果也好一点,实际上就是去掉dropout层,增加BN层,加速收敛控制过拟合等等的好处

    要注意由于用了数据增强的方法,输入就不再是RGB 三层的图像了,而是灰度图,所以预测的时候,读取图片需要这样做:
    img = load_img(img,grayscale=True)
    那么提取的血管总不能是灰色的吧,所以最后预测得到的(512,512)区间(0,1)的结果,乘回去原图就可以获得到血管了

    inputs = Input((self.img_rows, self.img_cols, 1))
    
    conv1 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(inputs)
    conv1 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    pool1 = layers.normalization.BatchNormalization()(pool1)
    
    conv2 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool1)
    conv2 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    pool2 = layers.normalization.BatchNormalization()(pool2)
    
    conv3 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool2)
    conv3 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    pool3 = layers.normalization.BatchNormalization()(pool3)
    
    conv4 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool3)
    conv4 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2))(conv4)
    pool4 = layers.normalization.BatchNormalization()(pool4)
    
    conv5 = Conv2D(1024, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool4)
    conv5 = Conv2D(1024, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv5)
    conv5 = layers.normalization.BatchNormalization()(conv5)
    
    up6 = Conv2D(512, 2, activation='relu', padding='same', kernel_initializer='he_normal')(
      UpSampling2D(size=(2, 2))(conv5))
    merge6 = concatenate([conv4, up6], axis=3)
    conv6 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge6)
    conv6 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv6)
    conv6 = layers.normalization.BatchNormalization()(conv6)
    
    up7 = Conv2D(256, 2, activation='relu', padding='same', kernel_initializer='he_normal')(
      UpSampling2D(size=(2, 2))(conv6))
    merge7 = concatenate([conv3, up7], axis=3)
    conv7 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge7)
    conv7 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv7)
    conv7 = layers.normalization.BatchNormalization()(conv7)
    
    up8 = Conv2D(128, 2, activation='relu', padding='same', kernel_initializer='he_normal')(
      UpSampling2D(size=(2, 2))(conv7))
    merge8 = concatenate([conv2, up8], axis=3)
    conv8 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge8)
    conv8 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv8)
    conv8 = layers.normalization.BatchNormalization()(conv8)
    
    up9 = Conv2D(64, 2, activation='relu', padding='same', kernel_initializer='he_normal')(
      UpSampling2D(size=(2, 2))(conv8))
    merge9 = concatenate([conv1, up9], axis=3)
    conv9 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge9)
    conv9 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv9)
    conv9 = layers.normalization.BatchNormalization()(conv9)
    
    conv9 = Conv2D(2, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv9)
    conv10 = Conv2D(1, 1, activation='sigmoid')(conv9)
    
    model = Model(inputs=inputs, outputs=conv10)
    model.compile(optimizer=Adam(lr=1e-4, loss="binary_crossentropy", metrics=['accuracy'])
    return model
    

    最后放一张实际预测的效果图,这个预测的准确度还是挺高的

    展开全文
    qq_31932741 2019-06-12 16:57:53
  • 67.97MB weixin_42166105 2021-05-01 20:17:08
  • 关于unet网络医学分割的网址 unet,大家可以在该网站中学习有关unet的知识 我将我的版本上传上了github,这是用keras实现的,运行data.py就可以将图片转换成.npy格式,这里我已经生成数据了,直接运行unet.py就...

    关于unet网络医学分割的网址 unet,大家可以在该网站中学习有关unet的知识

    我将我的版本上传上了github,这是用keras实现的,运行data.py就可以将图片转换成.npy格式,这里我已经生成数据了,直接运行unet.py就可以了

    大家可以在这里下载:https://github.com/huangshaoyin/U-net

    (具体参考的是哪个开源项目,我已经忘了,如有抄袭,可以联系我)

     

    打开data.py文件

    只要把这些路径改成你电脑对应的位置就行:

    现在已经改成了相对路径:

    还有unet.py

    这些路径也要修改,把路径都修改完了,就可以打开PyCharm运行unet.py

    训练图片:

    参考图片:

    测试图片:

    测试结果:

    测试结果很差,分析原因是因为训练样本只有30张,这是远远不够的,那篇论文说有做数据增强,以此来增加训练样本,所以接下来可以尝试一下增强数据后的训练结果。

     

    展开全文
    huangshaoyin 2018-07-14 10:32:05
  • 医学图像的数据格式十分复杂,数据形式有什么CT图像,MRI图像(目前还没接触过~),超声图像等(详细信息可以参考:医学图像了解。),而数据格式上有诸如dicom(.dcm,.IMA,常用的python处理库有:SimpleITK, ...

    这是一个学习记录博~可能有错,欢迎讨论

    P.S. 本文所用的unet源码来自Unet源码

    目标

    实现胃部超声图像的病灶分割

    医学数据以及预处理简介

    医学图像的数据格式十分复杂,数据形式有什么CT图像,MRI图像(目前还没接触过~),超声图像等(详细信息可以参考:医学图像了解。),而数据格式上有诸如dicom(.dcm,.IMA,常用的python处理库有:SimpleITK, pydicom库),Nifit(.nii或.nii.gz,常用的python处理库有:SimpleITK, nibabel库),NRRD(.nrrd,常用的python处理库有:SimpleITK,pynrrd库)等。

    这里再插一句,对于CT图像呢,有图像值与HU值这个概念,CT值属于医学领域的概念,通常称亨氏单位(hounsfield unit ,HU),反映了组织对X射线吸收程度,范围一般是[-1024,3071]。图像值属于计算机领域的概念,例如灰度图像的范围一般是[0,255]。

    在网上有看到这样的解释(CT图像之Hu值变换与窗宽窗位调整):

    对于nii格式的图片,经过测试,nibabel, simpleitk常用的api接口,都会自动的进行上述转化过程,即取出来的值已经是Hu了。(除非专门用nib.load('xx').dataobj.get_unscaled()或者itk.ReadImage('xx').GetPixel(x,y,z)才能取得原始数据)
    
    对于dcm格式的图片,经过测试,simpleitk, pydicom常用的api接口都不会将原始数据自动转化为Hu!!(itk snap软件读入dcm或nii都不会对数据进行scale操作)
    

    但我认为这个解释不是很对,因为我用simpleItk读取dcm格式的数据时,得出来的值就是Hu值(我也不太明白为什么),所以我觉得大家还是把数据打印出来看看,自己判断一下是不是需要转Hu值吧。

    对CT图像,在Hu值这个问题解决好之后,就要进行窗宽窗位的预处理了(窗宽窗位的解释参考来自CT图像之Hu值变换与窗宽窗位调整)。窗宽(window width)和窗位(window center),是用于选择感兴趣的CT值范围的,因为各种组织结构或病变具有不同的CT值,因此想要获得某一组织结构的最佳显示时,要选择适合观察该组织或病变的窗宽和窗位。

    举个例子,CT原图可能是这样的:
    在这里插入图片描述
    而在选取了窗宽窗位并对原图进行了clip之后的图是这样的:
    在这里插入图片描述
    是不是清晰了很多呢。

    本实验数据及预处理

    实验数据

    实验使用的数据为胃部超声图像,这里使用的超声数据是一般导出jpg格式的(一般会转成nrrd或者dcm来勾画病灶)。

    数据预处理

    增强对比度

    超声图像有的会比较暗,需要增加一下对比度的话可以使用如python+opencv直方图均衡化,以及tensorflow中常见的图片预处理操作等方法。

    数据增强

    数据量不够的情况下需要进行数据增强,可以使用keras中的Data generator来进行操作,具体做法可以参考Unet源码

    提取包含病灶的slices

    CT图像一般都很多的slices组成,很多时候病灶不会在全部的slices上出现,所以需要我们提取一下包含有病灶的slices,虽然本实验不涉及这个问题,但还是介绍一下吧。这个过程呢,就是检测一下mask数据,看一下在该层slice上是否有病灶出现,有就保留该slice,没有就舍弃。

    def getRangImageDepth(image):
        """
        args:
        image ndarray of shape (depth, height, weight)
        """
        # 得到轴向上出现过目标(label>=1)的切片
        ### image.shape为(depth,height,width),否则要改axis的值
        z = np.any(image, axis=(1,2)) # z.shape:(depth,)
        ###  np.where(z)[0] 输出满足条件,即true的元素位置,返回为一个列表,用0,-1取出第一个位置和最后一个位置的值
        startposition,endposition = np.where(z)[0][[0,-1]]
        return startposition, endposition
    

    Unet训练中需要注意的几个点

    这个Unet源码是用keras写的,代码比较简单易懂,所以这里不对代码做过多的介绍,如有问题,可以下面留言,这里只记录一下容易踩坑的点~

    图像的位深

    如果采用的数据是jpg格式的,可能会有图像位深问题的存在,从而导致预测失败,具体可以参考Unet训练自己的数据集,预测结果全黑

    检查数据类型

    在这个Unet源码的data.py文件中有对图像做一些预处理,比如reshape成网络所需的4维张量,以及对图像/255以将uint8转为网络所需的float型,这个根据自己的数据可能需要对代码进行一些修改,所以不要忘记检查这部分。

    输入图像大小的设置

    unet中会下采样4次,并且还要在上采样时还要与下采样中的特征图进行concatenate,因此,当图片的size在下采样时不能一直维持在偶数时,上采样进行特征结合就会出现问题。例如,37x37的下采样变成了18x18,而18x18在上采样时会变成36x36,在进行concatenate时,36x36和37x37就会出现大小不匹配的问题。

    对于这个问题,简单一点的解决方法就是把图片设置成每次下采样都能满足偶数的大小,如256x256;如果不想这样的话,我们也可以用tf.pad对图像张量进行填充,举例:

    ### 判断concatenate的特征图大小是否一致
    if drop4.get_shape().as_list() != up6.get_shape().as_list():
        # _,height, width, depth = up6.shape
        # 只要你使用Model,就必须保证该函数内全为layer,不能有其他函数,如果有其他函数必须用Lambda封装为layer
        up6_padded = Lambda(lambda x: tf.pad(x, [[0, 0], [1, 0], [0, 0], [0, 0]], 'REFLECT'))(up6)
        merge6 = concatenate([drop4, up6_padded], axis=3)
    else:
        merge6 = concatenate([drop4, up6], axis=3)
    

    这里对tf.pad进行一下简要说明,

    tf.pad(
        tensor,
        paddings,
        mode='CONSTANT',
        name=None
    )
    
    • tensor是待填充的张量
    • paddings指出要给tensor的哪个维度进行填充,以及填充方式,要注意的是paddings的rank必须和tensor的rank相同
    • mode指出用什么进行填充,可以取三个值,分别是"CONSTANT" ,“REFLECT”,“SYMMETRIC,mode=“CONSTANT” 填充0;mode="REFLECT"映射填充,按边缘翻且不复制边缘;mode="SYMMETRIC"对称填充,从对称轴开始复制,把边缘也复制了。
    • name是这个节点的名字

    在本实验中,我们的tensor是一个四维张量, [[0, 0], [1, 0], [0, 0], [0, 0]]分别对应[,height,width,depth],0表示不进行填充, [1, 0]表示对在图片上边填充一行,下边不填充,如果设置为[1,1]的话就代表在图片上下各填充1行;同样的,如果width为[2,3]就表示在图片的左边填充2列,右边填充3列。

    更换loss函数

    Unet源码中采用的是交叉熵损失来进行训练的,对于医学图像而言,有时候需要分割的病灶很小,这时候使用交叉熵损失有训练失败的可能,需要更换loss函数,这里有一个讲解的比较详细的博客,可以参考一下,从loss处理图像分割中类别极度不均衡的状况—keras

    例如将交叉熵损失换成Dice loss,具体做法为:

    def dice_coef(y_true, y_pred, smooth=1):
        intersection = K.sum(y_true * y_pred, axis=[1,2,3])
        union = K.sum(y_true, axis=[1,2,3]) + K.sum(y_pred, axis=[1,2,3])
        return K.mean( (2. * intersection + smooth) / (union + smooth), axis=0)
    
    
    def dice_coef_loss(y_true, y_pred):
    
        return 1 - dice_coef(y_true, y_pred, smooth=1)
    

    然后将modle.compile改一下:

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

    更改学习率、batchsize和epoch

    学习率不对也可能会导致预测失败,所以需要自己调整,然后在测试看看。
    还有batchsize和epoch也要根据自己的数据集的情况进行调整。

    emmm…暂时就这么多啦~

    展开全文
    sinat_35779431 2020-05-12 09:39:51
  • 这次实验用来学习unet网络实现图像分割(keras, backend: tensorflow)。 数据集DRIVE:为眼部图像,目的是分割出眼部血管。 数据集结构: 上面分别是训练的原始图片images、first_manual、mask 整体流程: 1、前期...

    这次实验用来学习unet网络实现图像分割(keras, backend: tensorflow)。
    数据集DRIVE:为眼部图像,目的是分割出眼部血管。
    数据集结构:在这里插入图片描述
    train_imagesfirst_manual
    在这里插入图片描述
    上面分别是训练的原始图片images、first_manual、mask

    整体流程:

    • 1、前期准备:将所有图片写入h5文件,(为了后期训练网络时减少io时间)
    • 2、训练网络
      • 2.1、获取训练图片和groundtruth图片
        • 2.1.1、 读取hdf5文件,获取img array
        • 2.1.2、 预处理所有的img
        • 2.2.3、 提取patches
      • 2.2、构建model,并训练
      • 2.3、模型保存(包括architecture 和 weights)
    • 3、测试
      • 3.1、获取测试图片
      • 3.2、model加载
      • 3.3、model predict
      • 3.4、测试结果重建,可视化显示分割结果(patches->full image)
      • 3.5、用不同的metrics测评model结果

    文档结构
    在这里插入图片描述
    注意,lib文件夹里的加上__init__.py很重要,加上后,这个文件夹才能被python识别为一个package,从而,可以在run_training.py中写from lib.help_function import *

    第一步
    python prepare_dataset_DRIVE.py
    将图片读取为ndarray格式,并保存至hdf5文件

    # -*- coding: utf-8 -*-
    
    #==========================================================
    #
    #  This prepare the hdf5 datasets of the DRIVE database
    #
    #============================================================
    
    import os
    import h5py
    import numpy as np
    from PIL import Image
    """
    PIL 读取的图片为RGB通道的Image格式,
    cv2 读取的图片为BGR通道的ndarray格式
    skimage 读取的图片为RGB通道的ndarray格式
    """
    
    
    
    def write_hdf5(arr,outfile):
      with h5py.File(outfile,"w") as f:
        f.create_dataset("image", data=arr, dtype=arr.dtype)
    
    
    #------------Path of the images --------------------------------------------------------------
    #train
    original_imgs_train = "./DRIVE/training/images/"
    groundTruth_imgs_train = "./DRIVE/training/1st_manual/"
    borderMasks_imgs_train = "./DRIVE/training/mask/"
    #test
    original_imgs_test = "./DRIVE/test/images/"
    groundTruth_imgs_test = "./DRIVE/test/1st_manual/"
    borderMasks_imgs_test = "./DRIVE/test/mask/"
    #---------------------------------------------------------------------------------------------
    
    # 图像原始数量、通道数、高、宽
    Nimgs = 20
    channels = 3
    height = 584
    width = 565
    dataset_path = "./DRIVE_datasets_training_testing/"
    
    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)")
       
        assert(imgs.shape == (Nimgs,height,width,channels))
        groundTruth = np.reshape(groundTruth,(Nimgs,height,width,1))
        border_masks = np.reshape(border_masks,(Nimgs,height,width,1))
        assert(groundTruth.shape == (Nimgs,height,width,1))
        assert(border_masks.shape == (Nimgs,height,width,1))
        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")
    

    配置文件configuration.txt

    [data paths]
    path_local =  ./DRIVE_datasets_training_testing/
    train_imgs_original = DRIVE_dataset_imgs_train.hdf5
    train_groundTruth = DRIVE_dataset_groundTruth_train.hdf5
    train_border_masks = DRIVE_dataset_borderMasks_train.hdf5
    test_imgs_original = DRIVE_dataset_imgs_test.hdf5
    test_groundTruth = DRIVE_dataset_groundTruth_test.hdf5
    test_border_masks = DRIVE_dataset_borderMasks_test.hdf5
    
    
    
    [experiment name]
    name = test
    
    
    [data attributes]
    #Dimensions of the patches extracted from the full images
    patch_height = 48
    patch_width = 48
    
    
    [training settings]
    #number of total patches:
    N_subimgs = 190000
    #if patches are extracted only inside the field of view:
    inside_FOV = False
    #Number of training epochs
    N_epochs = 150
    batch_size = 32
    #if running with nohup
    nohup = True
    
    
    [testing settings]
    #Choose the model to test: best==epoch with min loss, last==last epoch
    best_last = best
    #number of full images for the test (max 20)
    full_images_to_test = 20
    #How many original-groundTruth-prediction images are visualized in each image
    N_group_visual = 1
    #Compute average in the prediction, improve results but require more patches to be predicted
    average_mode = True
    #Only if average_mode==True. Stride for patch extraction, lower value require more patches to be predicted
    stride_height = 5
    stride_width = 5
    #if running with nohup
    nohup = False
    

    第二步
    python run_training.py

    # -*- coding: utf-8 -*-
    
    ###################################################
    #
    #   run_training.py  Script to launch the training
    #
    ##################################################
    
    import os, sys
    # 注意,py3是configparser,py2是ConfigParser
    import configparser
    
    
    #config file to read from
    config = configparser.RawConfigParser()
    config.readfp(open(r'./configuration.txt'))
    #===========================================
    #name of the experiment
    name_experiment = config.get('experiment name', 'name')
    nohup = config.getboolean('training settings', 'nohup')   #std output on log file?
    
    # 在自己机器上测试时,注释掉
    run_GPU = '' if sys.platform == 'win32' else ' THEANO_FLAGS=device=gpu,floatX=float32 '
    
    #create a folder for the results
    result_dir = name_experiment
    print("\n1. Create directory for the results (if not already existing)")
    if os.path.exists(result_dir):
        print("Dir already existing")
    elif sys.platform=='win32':
        os.system('mkdir ' + result_dir)
    else:
        os.system('mkdir -p ' +result_dir)
    
    print("copy the configuration file in the results folder")
    if sys.platform=='win32':
        os.system('copy configuration.txt .\\' +name_experiment+'\\'+name_experiment+'_configuration.txt')
    else:
        os.system('cp configuration.txt ./' +name_experiment+'/'+name_experiment+'_configuration.txt')
    
    # 注意,在自己机器上测试时,把与run_GPU有关的命令去掉
    # run the experiment
    if nohup:
        print("\n2. Run the training on GPU with nohup")
        #os.system(run_GPU +' nohup python -u ./src/retinaNN_training.py > ' +'./'+name_experiment+'/'+name_experiment+'_training.nohup')
        os.system(run_GPU + 'nohup python -u ./src/retinaNN_training.py > ' +'./'+name_experiment+'/'+name_experiment+'_training.nohup')
    
    else:
        print("\n2. Run the training on GPU (no nohup)")
        os.system(run_GPU + 'python ./src/retinaNN_training.py')
    
    # os.system('python ./src/retinaNN_training.py')
    
    # -*- coding: utf-8 -*-
    # retinaNN_training.py
    import numpy as np
    import configparser
    
    from keras.models import Model
    from keras.layers import Input, concatenate, Conv2D, MaxPooling2D, UpSampling2D, Reshape, core, Dropout
    from keras.optimizers import Adam
    from keras.callbacks import ModelCheckpoint, LearningRateScheduler
    from keras import backend as K
    from keras.utils.vis_utils import plot_model as plot
    from keras.optimizers import SGD
    
    import sys
    sys.path.insert(0, './lib/') 
    from help_functions import *
    
    #function to obtain data for training/testing (validation)
    from extract_patches import get_data_training
    
    
    
    #Define the neural network
    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)
        #  注意,写成这种结构,并且用的loss为categorical_crossentropy,
        # 需要对groundtruth数据进行处理,见后面help_function.py里的mask_Unet
        conv6 = Conv2D(2, (1, 1), activation='relu',padding='same')(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
    
    
    
    
    
    print("read configure")
    #========= Load settings from Config file
    config = configparser.RawConfigParser()
    config.read('configuration.txt')
    #patch to the datasets
    path_data = config.get('data paths', 'path_local')
    #Experiment name
    name_experiment = config.get('experiment name', 'name')
    #training settings
    N_epochs = int(config.get('training settings', 'N_epochs'))
    batch_size = int(config.get('training settings', 'batch_size'))
    
    
    print("load data")
    #============ Load the data and divided in patches
    patches_imgs_train, patches_masks_train = get_data_training(
        DRIVE_train_imgs_original = path_data + config.get('data paths', 'train_imgs_original'),
        DRIVE_train_groudTruth = path_data + config.get('data paths', 'train_groundTruth'),  #masks
        patch_height = int(config.get('data attributes', 'patch_height')),
        patch_width = int(config.get('data attributes', 'patch_width')),
        N_subimgs = int(config.get('training settings', 'N_subimgs')),
        inside_FOV = config.getboolean('training settings', 'inside_FOV') #select the patches only inside the FOV  (default == True)
    )
    
    # 这些可以不写
    print("create sampel")
    #========= Save a sample of what you're feeding to the neural network ==========
    N_sample = min(patches_imgs_train.shape[0],40)
    visualize(group_images(patches_imgs_train[0:N_sample,:,:,:],5),'./'+name_experiment+'/'+"sample_input_imgs")#.show()
    visualize(group_images(patches_masks_train[0:N_sample,:,:,:],5),'./'+name_experiment+'/'+"sample_input_masks")#.show()
    
    print("construct model")
    #=========== Construct and save the model arcitecture =====
    n_ch = patches_imgs_train.shape[3]
    patch_height = patches_imgs_train.shape[1]
    patch_width = patches_imgs_train.shape[2]
    model = get_unet(n_ch, patch_height, patch_width)  #the U-net model
    print("Check: final output of the network:")
    print(model.output_shape)
    plot(model, to_file='./'+name_experiment+'/'+name_experiment + '_model.png')   #check how the model looks like
    json_string = model.to_json()
    open('./'+name_experiment+'/'+name_experiment +'_architecture.json', 'w').write(json_string)
    
    #============  Training ==================================
    checkpointer = ModelCheckpoint(filepath='./'+name_experiment+'/'+name_experiment +'_best_weights.h5', verbose=1, monitor='val_loss', mode='auto', save_best_only=True) #save at each epoch if the validation decreased
    
    
    # def step_decay(epoch):
    #     lrate = 0.01 #the initial learning rate (by default in keras)
    #     if epoch==100:
    #         return 0.005
    #     else:
    #         return lrate
    #
    # lrate_drop = LearningRateScheduler(step_decay)
    
    patches_masks_train = masks_Unet(patches_masks_train)  #reduce memory consumption
    model.fit(patches_imgs_train, patches_masks_train, nb_epoch=N_epochs, batch_size=batch_size, verbose=1, shuffle=True, validation_split=0.1, callbacks=[checkpointer])
    
    
    #========== Save and test the last model ===================
    model.save_weights('./'+name_experiment+'/'+name_experiment +'_last_weights.h5', overwrite=True)
    

    上面出现过的一些函数:

    # extract_patches.py
    import numpy as np
    import random
    import configparser
    
    from help_functions import load_hdf5
    from help_functions import visualize
    from help_functions import group_images
    
    from pre_processing import my_PreProc
    
    
    #To select the same images
    # random.seed(10)
    
    #Load the original data and return the extracted patches for training/testing
    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)
        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
    
    
    #data consinstency check
    def data_consistency_check(imgs,masks):
        assert(len(imgs.shape)==len(masks.shape))
        assert(imgs.shape[0]==masks.shape[0])
        assert(imgs.shape[1]==masks.shape[1])
        assert(imgs.shape[2]==masks.shape[2])
        assert(masks.shape[3]==1)
        assert(imgs.shape[3]==1 or imgs.shape[3]==3)
    
    
    #extract patches randomly in the full training images
    #  -- Inside OR in full image
    def extract_random(full_imgs,full_masks, patch_h,patch_w, N_patches, inside=True):
        if (N_patches%full_imgs.shape[0] != 0):
            print("N_patches: plase enter a multiple of 20")
            exit()
        assert (len(full_imgs.shape)==4 and len(full_masks.shape)==4)  #4D arrays
        assert (full_imgs.shape[3]==1 or full_imgs.shape[3]==3)  #check the channel is 1 or 3
        assert (full_masks.shape[3]==1)   #masks only black and white
        assert (full_imgs.shape[2] == full_masks.shape[2] and full_imgs.shape[1] == full_masks.shape[1])
        patches = np.empty((N_patches,patch_h,patch_w,full_imgs.shape[3]))
        patches_masks = np.empty((N_patches,patch_h,patch_w,full_masks.shape[3]))
        img_h = full_imgs.shape[1]  #height of the full image
        img_w = full_imgs.shape[2] #width of the full image
        # (0,0) in the center of the image
        patch_per_img = int(N_patches/full_imgs.shape[0])  #N_patches equally divided in the full images
        print("patches per full image: ", str(patch_per_img))
        iter_tot = 0   #iter over the total numbe rof patches (N_patches)
        for i in range(full_imgs.shape[0]):  #loop over the full images
            k=0
            while k <patch_per_img:
                x_center = random.randint(0+int(patch_w/2),img_w-int(patch_w/2))
                # print "x_center " +str(x_center)
                y_center = random.randint(0+int(patch_h/2),img_h-int(patch_h/2))
                # print "y_center " +str(y_center)
                #check whether the patch is fully contained in the FOV
                if inside==True:
                    if is_patch_inside_FOV(x_center,y_center,img_w,img_h,patch_h)==False:
                        continue
                patch = full_imgs[i,y_center-int(patch_h/2):y_center+int(patch_h/2),x_center-int(patch_w/2):x_center+int(patch_w/2),:]
                patch_mask = full_masks[i,y_center-int(patch_h/2):y_center+int(patch_h/2),x_center-int(patch_w/2):x_center+int(patch_w/2),:]
                patches[iter_tot]=patch
                patches_masks[iter_tot]=patch_mask
                iter_tot +=1   #total
                k+=1  #per full_img
        return patches, patches_masks
    
    #check if the patch is fully contained in the FOV
    def is_patch_inside_FOV(x,y,img_w,img_h,patch_h):
        x_ = x - int(img_w/2) # origin (0,0) shifted to image center
        y_ = y - int(img_h/2)  # origin (0,0) shifted to image center
        R_inside = 270 - int(patch_h * np.sqrt(2.0) / 2.0) #radius is 270 (from DRIVE db docs), minus the patch diagonal (assumed it is a square #this is the limit to contain the full patch in the FOV
        radius = np.sqrt((x_*x_)+(y_*y_))
        if radius < R_inside:
            return True
        else:
            return False
    

    help_function.py

    # -*- coding: utf-8 -*-
    
    import h5py
    import numpy as np
    from PIL import Image
    from matplotlib import pyplot as plt
    
    def load_hdf5(infile):
      with h5py.File(infile,"r") as f:  #"with" close the file after its nested commands
        return f["image"][()]
    
    def write_hdf5(arr,outfile):
      with h5py.File(outfile,"w") as f:
        f.create_dataset("image", data=arr, dtype=arr.dtype)
    
    #convert RGB image in black and white
    def rgb2gray(rgb):
        assert (len(rgb.shape)==4)  #4D arrays
        assert (rgb.shape[3]==3)
        bn_imgs = rgb[:,:,:,0]*0.299 + rgb[:,:,:,1]*0.587 + rgb[:,:,:,2]*0.114
        bn_imgs = np.reshape(bn_imgs,(rgb.shape[0],rgb.shape[1],rgb.shape[2],1))
        return bn_imgs
    
    #group a set of images row per columns
    def group_images(data,per_row):
        assert data.shape[0]%per_row==0
        assert (data.shape[3]==1 or data.shape[3]==3)
        all_stripe = []
        for i in range(int(data.shape[0]/per_row)):
            stripe = data[i*per_row]
            for k in range(i*per_row+1, i*per_row+per_row):
                stripe = np.concatenate((stripe,data[k]),axis=1)
            all_stripe.append(stripe)
        totimg = all_stripe[0]
        for i in range(1,len(all_stripe)):
            totimg = np.concatenate((totimg,all_stripe[i]),axis=0)
        return totimg
    
    
    #visualize image (as PIL image, NOT as matplotlib!)
    def visualize(data,filename):
        assert (len(data.shape)==3) #height*width*channels
        img = None
        if data.shape[2]==1:  #in case it is black and white
            data = np.reshape(data,(data.shape[0],data.shape[1]))
        if np.max(data)>1:
            img = Image.fromarray(data.astype(np.uint8))   #the image is already 0-255
        else:
            img = Image.fromarray((data*255).astype(np.uint8))  #the image is between 0-1
        img.save(filename + '.png')
        return img
    
    
    #prepare the mask in the right shape for the Unet
    def masks_Unet(masks):
        assert (len(masks.shape)==4)  #4D arrays
        assert (masks.shape[3]==1 )  #check the channel is 1
        im_h = masks.shape[1]
        im_w = masks.shape[2]
        masks = np.reshape(masks,(masks.shape[0],im_h*im_w))
        new_masks = np.empty((masks.shape[0],im_h*im_w,2))
        for i in range(masks.shape[0]):
            for j in range(im_h*im_w):
                if  masks[i,j] == 0:
                    new_masks[i,j,0]=1
                    new_masks[i,j,1]=0
                else:
                    new_masks[i,j,0]=0
                    new_masks[i,j,1]=1
        return new_masks
    
    
    def pred_to_imgs(pred, patch_height, patch_width, mode="original"):
        assert (len(pred.shape)==3)  #3D array: (Npatches,height*width,2)
        assert (pred.shape[2]==2 )  #check the classes are 2
        pred_images = np.empty((pred.shape[0],pred.shape[1]))  #(Npatches,height*width)
        if mode=="original":
            for i in range(pred.shape[0]):
                for pix in range(pred.shape[1]):
                    pred_images[i,pix]=pred[i,pix,1]
        elif mode=="threshold":
            for i in range(pred.shape[0]):
                for pix in range(pred.shape[1]):
                    if pred[i,pix,1]>=0.5:
                        pred_images[i,pix]=1
                    else:
                        pred_images[i,pix]=0
        else:
            print("mode ", str(mode), " not recognized, it can be 'original' or 'threshold'")
            exit()
        pred_images = np.reshape(pred_images,(pred_images.shape[0], patch_height, patch_width,1))
        return pred_images
    

    pre_processing.py

    # -*- coding: utf-8 -*-
    
    ###################################################
    #
    #  pre_processing.py Script to pre-process the original imgs
    #
    ##################################################
    
    
    import numpy as np
    from PIL import Image
    import cv2
    
    from help_functions import *
    
    
    #My pre processing (use for both training and testing!)
    def my_PreProc(data):
        assert(len(data.shape)==4)
        assert (data.shape[3]==3)  #Use the original images
        #black-white conversion
        train_imgs = rgb2gray(data)
        #my preprocessing:
        train_imgs = dataset_normalized(train_imgs)
        train_imgs = clahe_equalized(train_imgs)
        train_imgs = adjust_gamma(train_imgs, 1.2)
        train_imgs = train_imgs/255.  #reduce to 0-1 range
        return train_imgs
    
    
    #============================================================
    #========= PRE PROCESSING FUNCTIONS ========================#
    #============================================================
    
    #==== histogram equalization
    def histo_equalized(imgs):
        assert (len(imgs.shape)==4)  #4D arrays
        assert (imgs.shape[3]==1)  #check the channel is 1
        imgs_equalized = np.empty(imgs.shape)
        for i in range(imgs.shape[0]):
            imgs_equalized[i,:,:,0] = cv2.equalizeHist(np.array(imgs[i,:,:,0], dtype = np.uint8))
        return imgs_equalized
    
    
    # CLAHE (Contrast Limited Adaptive Histogram Equalization)
    #adaptive histogram equalization is used. In this, image is divided into small blocks called "tiles" (tileSize is 8x8 by default in OpenCV). Then each of these blocks are histogram equalized as usual. So in a small area, histogram would confine to a small region (unless there is noise). If noise is there, it will be amplified. To avoid this, contrast limiting is applied. If any histogram bin is above the specified contrast limit (by default 40 in OpenCV), those pixels are clipped and distributed uniformly to other bins before applying histogram equalization. After equalization, to remove artifacts in tile borders, bilinear interpolation is applied
    def clahe_equalized(imgs):
        assert (len(imgs.shape)==4)  #4D arrays
        assert (imgs.shape[3]==1)  #check the channel is 1
        #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(imgs.shape[0]):
            imgs_equalized[i,:,:,0] = clahe.apply(np.array(imgs[i,:,:,0], dtype = np.uint8))
        return imgs_equalized
    
    
    # ===== normalize over the dataset
    def dataset_normalized(imgs):
        assert (len(imgs.shape)==4)  #4D arrays
        assert (imgs.shape[3]==1)  #check the channel is 1
        imgs_normalized = np.empty(imgs.shape)
        imgs_std = np.std(imgs)
        imgs_mean = np.mean(imgs)
        imgs_normalized = (imgs-imgs_mean)/imgs_std
        for i in range(imgs.shape[0]):
            imgs_normalized[i] = ((imgs_normalized[i] - np.min(imgs_normalized[i])) / (np.max(imgs_normalized[i])-np.min(imgs_normalized[i])))*255
        return imgs_normalized
    
    
    def adjust_gamma(imgs, gamma=1.0):
        assert (len(imgs.shape)==4)  #4D arrays
        assert (imgs.shape[3]==1)  #check the channel is 1
        # build a lookup table mapping the pixel values [0, 255] to
        # their adjusted gamma values
        invGamma = 1.0 / gamma
        table = np.array([((i / 255.0) ** invGamma) * 255 for i in np.arange(0, 256)]).astype("uint8")
        # apply gamma correction using the lookup table
        new_imgs = np.empty(imgs.shape)
        for i in range(imgs.shape[0]):
            new_imgs[i,:,:,0] = cv2.LUT(np.array(imgs[i,:,:,0], dtype = np.uint8), table)
        return new_imgs
    
    def pred_to_imgs(pred, patch_height, patch_width, mode="original"):
        assert (len(pred.shape)==3)  #3D array: (Npatches,height*width,2)
        assert (pred.shape[2]==2 )  #check the classes are 2
        pred_images = np.empty((pred.shape[0],pred.shape[1]))  #(Npatches,height*width)
        if mode=="original":
            for i in range(pred.shape[0]):
                for pix in range(pred.shape[1]):
                    pred_images[i,pix]=pred[i,pix,1]
        elif mode=="threshold":
            for i in range(pred.shape[0]):
                for pix in range(pred.shape[1]):
                    if pred[i,pix,1]>=0.5:
                        pred_images[i,pix]=1
                    else:
                        pred_images[i,pix]=0
        else:
            print("mode ", str(mode), " not recognized, it can be 'original' or 'threshold'")
            exit()
        pred_images = np.reshape(pred_images,(pred_images.shape[0], patch_height, patch_width, 1))
        return pred_images
    

    训练的步数少时,出来的效果如下,我还以为是网络写错了:
    在这里插入图片描述
    随着训练步数再增加一点,慢慢有了雏形:
    在这里插入图片描述
    训练步数再增加一点,可见继续训练下去,是能达到一定的效果的,这次只是为了学习,就只跑了一小会:
    在这里插入图片描述

    展开全文
    normol 2019-02-27 17:39:35
  • weixin_44485421 2021-08-16 10:50:09
  • jizhidexiaoming 2019-07-23 11:03:57
  • weixin_47035018 2021-06-19 10:29:07
  • 60KB weixin_42144201 2021-05-04 18:43:02
  • qq_35793394 2020-06-23 13:48:56
  • hehuaiyuyu 2020-05-23 23:06:49
  • weixin_45435630 2020-06-18 22:37:56
  • flyfor2013 2019-10-10 10:00:00
  • shajiayu1 2020-12-20 22:13:42
  • qq_42722197 2021-06-12 00:56:07
  • weixin_42097914 2021-03-07 06:03:36
  • qq_41375318 2020-04-10 21:26:23
  • 160KB weixin_42144201 2021-05-11 13:17:44
  • 5.66MB weixin_42097369 2021-06-18 19:13:18
  • 16.37MB weixin_42159267 2021-02-14 05:58:39
  • xholes 2019-05-10 07:44:18
  • normol 2019-03-24 22:14:55
  • guzhao9901 2021-08-09 18:13:45
  • algorithmPro 2020-09-26 11:40:00
  • 3星
    47.26MB weixin_39841365 2019-08-11 06:35:34
  • 468KB k3v1n1990s 2018-01-04 00:36:05

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,329
精华内容 531
关键字:

unet医学图像分割