Unet实现图像分割（四）

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 = 72
patch_width = 72[training settings]
#number of total patches:
N_subimgs = 38000
#if patches are extracted only inside the field of view:
inside_FOV = False
#Number of training epochs
N_epochs = 20
batch_size = 32
#if running with nohup
nohup = False[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 = 20
stride_width = 20
#if running with nohup
nohup = False

prepare_datasets_DRIVE.py

#==========================================================
#
#  This prepare the hdf5 datasets of the DRIVE database
#
#============================================================
import os
import h5py
import numpy as np
from PIL import Imagedef write_hdf5(arr,outfile): # 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))  # 二值图像 channels=1border_masks = np.empty((Nimgs,height,width)) # 二值图像 channels=1for path, subdirs, files in os.walk(imgs_dir): #list all files, directories in the path # path=当前路径 subdirs=子文件夹 files=文件夹内所有的文件for i in range(len(files)):  # len(files) 所有图像的数量#originalprint("original image: " +files[i])img = Image.open(imgs_dir+files[i]) # 读取图像到内存imgs[i] = np.asarray(img)           # 转换成numpy数据格式#corresponding ground truthgroundTruth_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 masksborder_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# 调整张量格式 [Nimg, channels, height, width]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_masksif 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")

retinaNN_training.py

###################################################
#
#   Script to:
#   - Load the images and extract the patches
#   - Define the neural network
#   - define the training
#
##################################################
import numpy as np
import configparser # Python 3.6中 configparser全使用小写
import os # os模块中主要用于处理文件和目录
from keras.models import Model
from keras.layers import Input, concatenate, Conv2D, MaxPooling2D, UpSampling2D, Reshape, core, Dropout #core内部定义了一系列常用的网络层，包括全连接、激活层等
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 * # 导入help_functions脚本文件中的所有函数#function to obtain data for training/testing (validation)
from extract_patches import get_data_training # 导入extract_patches 脚本中的 get_data_training函数#Define the neural network
def get_unet(n_ch,patch_height,patch_width):inputs = Input(shape=(n_ch,patch_height,patch_width))#data_format：字符串，“channels_first”或“channels_last”之一，代表图像的通道维的位置。#以128x128的RGB图像为例，“channels_first”应将数据组织为（3,128,128），而“channels_last”应将数据组织为（128,128,3）。该参数的默认值是~/.keras/keras.json中设置的值，若从未设置过，则为“channels_last”。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)##1×1的卷积的作用#大概有两个方面的作用：1. 实现跨通道的交互和信息整合2. 进行卷积核通道数的降维和升维。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=Adam(lr=0.001), loss='categorical_crossentropy',metrics=['accuracy'])return model'''模型Model的compile方法:compile(self, optimizer, loss, metrics=None, loss_weights=None, sample_weight_mode=None, weighted_metrics = None, target_tensors=None)本函数编译模型以供训练，参数有optimizer：         优化器，为预定义优化器名或优化器对.可以在调用model.compile()之前初始化一个优化器对象，然后传入该函数。loss：              损失函数，为预定义损失函数名或一个目标函数metrics：           列表，包含评估模型在训练和测试时的性能的指标，典型用法是metrics=['accuracy']如果要在多输出模型中为不同的输出指定不同的指标，可像该参数传递一个字典，例如metrics={'ouput_a': 'accuracy'}sample_weight_mode：如果需要按时间步为样本赋权（2D权矩阵），将该值设为“temporal”。默认为“None”，代表按样本赋权（1D权）。如果模型有多个输出，可以向该参数传入指定sample_weight_mode的字典或列表。在下面fit函数的解释中有相关的参考内容。weighted_metrics:   metrics列表，在训练和测试过程中，这些metrics将由sample_weight或clss_weight计算并赋权target_tensors:     默认情况下，Keras将为模型的目标创建一个占位符，该占位符在训练过程中将被目标数据代替。如果你想使用自己的目标张量（相应的，Keras将不会在训练时期望为这些目标张量载入外部的numpy数据），你可以通过该参数手动指定。目标张量可以是一个单独的张量（对应于单输出模型），也可以是一个张量列表，或者一个name->tensor的张量字典。kwargs：            使用TensorFlow作为后端请忽略该参数，若使用Theano/CNTK作为后端，kwargs的值将会传递给 K.function。如果使用TensorFlow为后端，这里的值会被传给tf.Session.run在Keras中，compile主要完成损失函数和优化器的一些配置，是为训练服务的。'''
#========= 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')) #训练的批量大小#============ 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'),  #maskspatch_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)
)
#========= 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()#=========== Construct and save the model arcitecture =====
#调用网络 及 保存网络模型
n_ch = patches_imgs_train.shape[1]
patch_height = patches_imgs_train.shape[2]
patch_width = patches_imgs_train.shape[3]
#U-net 网络 [batchsize, channels, patch_heigh, patch_width]
model = get_unet(n_ch, patch_height, patch_width)  #the U-net model
print("Check: final output of the network:")
print(model.output_shape)
#os.environ["PATH"] += os.pathsep + 'C:/Program Files (x86)/Graphviz2.38/bin/'
#调用pydot显示模型
#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 ==================================
#采用回调函数的形式保存每个epoch数据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
'''
keras.callbacks.ModelCheckpoint(filepath,monitor='val_loss', verbose=0, save_best_only=False, save_weights_only = False, mode='auto', period=1)
该回调函数将在每个epoch后保存模型到filepath;
filename：字符串，保存模型的路径
monitor：需要监视的值
verbose：信息展示模式，0或1
save_best_only：当设置为True时，将只保存在验证集上性能最好的模型
mode：‘auto’，‘min’，‘max’之一，在save_best_only=True时决定性能最佳模型的评判准则，例如，当监测值为val_acc时，模式应为max，当检测值为val_loss时，模式应为min。在auto模式下，评价准则由被监测值的名字自动推断。
save_weights_only：若设置为True，则只保存模型权重，否则将保存整个模型（包括模型结构，配置信息等）
period：CheckPoint之间的间隔的epoch数
'''# 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)
'''
keras.callbacks.LearningRateScheduler(schedule)
该回调函数是学习率调度器
schedule：函数，该函数以epoch号为参数（从0算起的整数），返回一个新学习率（浮点数）
'''patches_masks_train = masks_Unet(patches_masks_train)  # reduce memory consumption
model.fit(patches_imgs_train, patches_masks_train, epochs=N_epochs, batch_size=batch_size, verbose=1, shuffle=True, validation_split=0.1, callbacks=[checkpointer])'''
model.fit(self, x, y, batch_size=32, epochs=10, verbose=1, callbacks=None, validation_split=0.0, validation_data=None, shuffle=True, class_weight=None, sample_weight=None, initial_epoch=0)
本函数将模型训练nb_epoch轮，其参数有：
x：输入数据。如果模型只有一个输入，那么x的类型是numpy array，如果模型有多个输入，那么x的类型应当为list，list的元素是对应于各个输入的numpy array
y：标签，numpy array
batch_size：整数，指定进行梯度下降时每个batch包含的样本数。
epochs：整数，训练终止时的epoch值，训练将在达到该epoch值时停止
verbose：日志显示，0为不在标准输出流输出日志信息，1为输出进度条记录，2为每个epoch输出一行记录
callbacks：list，其中的元素是keras.callbacks.Callback的对象，如ModelCheckpoint()、LearningRateScheduler()等。
validation_split：0~1之间的浮点数，用来指定训练集的一定比例数据作为验证集。注意，validation_split的划分在shuffle之前，因此如果你的数据本身是有序的，需要先手工打乱再指定validation_split，否则可能会出现验证集样本不均匀。
validation_data：形式为（X，y）的tuple，是指定的验证集。此参数将覆盖validation_spilt。
shuffle：布尔值或字符串，一般为布尔值，表示是否在训练过程中随机打乱输入样本的顺序。若为字符串“batch”，则是用来处理HDF5数据的特殊情况，它将在batch内部将数据打乱。
class_weight：字典，将不同的类别映射为不同的权值，该参数用来在训练过程中调整损失函数（只能用于训练）
sample_weight：权值的numpy array，用于在训练时调整损失函数（仅用于训练）。可以传递一个1D的与样本等长的向量用于对样本进行1对1的加权，或者在面对时序数据时，传递一个的形式为（samples，sequence_length）的矩阵来为每个时间步上的样本赋不同的权。这种情况下请确定在编译模型时添加了sample_weight_mode='temporal'。
initial_epoch: 从该参数指定的epoch开始训练，在继续之前的训练时有用。
'''#========== Save and test the last model ===================
model.save_weights('./'+name_experiment+'/'+name_experiment +'_last_weights.h5', overwrite=True)
'''
Keras中模型的保存分为两部分分别是保存架构jasonfile.write()和权重save_weights();同时模型的读取也包括网络架构读取model = model_from_json(open('').read()) 和模型训练好的权重读取model.load_weights('')。
'''
#test the model
# score = model.evaluate(patches_imgs_test, masks_Unet(patches_masks_test), verbose=0)
# print('Test score:', score[0])
# print('Test accuracy:', score[1])

retinaNN_predict.py

###################################################
#
#   Script to
#   - Calculate prediction of the test dataset
#   - Calculate the parameters to evaluate the prediction
#
##################################################
#Python
import numpy as np
import configparser
from matplotlib import pyplot as plt
#Keras
from keras.models import model_from_json
from keras.models import Model
#scikit learn
#sklearn.metric提供了一些函数，用来计算真实值与预测值之间的预测误差；这里用的评价标准主要集中如下几个方面：
from sklearn.metrics import roc_curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import jaccard_similarity_score
from sklearn.metrics import f1_score
import sys
sys.path.insert(0, './lib/')
# help_functions.py
from help_functions import *
# extract_patches.py
from extract_patches import recompone
from extract_patches import recompone_overlap
from extract_patches import paint_border
from extract_patches import kill_border
from extract_patches import pred_only_FOV
from extract_patches import get_data_testing
from extract_patches import get_data_testing_overlap
# pre_processing.py
from pre_processing import my_PreProc#========= CONFIG FILE TO READ FROM =======
config = configparser.RawConfigParser()
config.read('configuration.txt')
#===========================================
#run the training on invariant or local
path_data = config.get('data paths', 'path_local')  #数据路径#original test images (for FOV selection)
DRIVE_test_imgs_original = path_data + config.get('data paths', 'test_imgs_original')  #测试集图像封装文件
test_imgs_orig = load_hdf5(DRIVE_test_imgs_original) #测试集图像
full_img_height = test_imgs_orig.shape[2]
full_img_width = test_imgs_orig.shape[3]
#the border masks provided by the DRIVE
DRIVE_test_border_masks = path_data + config.get('data paths', 'test_border_masks')  #测试集掩膜封装文件
test_border_masks = load_hdf5(DRIVE_test_border_masks)
# dimension of the patches
# 图像块的维度
patch_height = int(config.get('data attributes', 'patch_height'))
patch_width = int(config.get('data attributes', 'patch_width'))
#the stride in case output with average
# 图像分块的跳跃步长
stride_height = int(config.get('testing settings', 'stride_height'))
stride_width = int(config.get('testing settings', 'stride_width'))
assert (stride_height < patch_height and stride_width < patch_width)
#model name
name_experiment = config.get('experiment name', 'name')
path_experiment = './' +name_experiment +'/'
#N full images to be predicted
Imgs_to_test = int(config.get('testing settings', 'full_images_to_test')) # 20张图像全部进行预测
#Grouping of the predicted images
N_visual = int(config.get('testing settings', 'N_group_visual'))
#====== average mode ===========
average_mode = config.getboolean('testing settings', 'average_mode')  # average=True# #ground truth
# gtruth= path_data + config.get('data paths', 'test_groundTruth')
# img_truth= load_hdf5(gtruth)
# visualize(group_images(test_imgs_orig[0:20,:,:,:],5),'original')#.show()
# visualize(group_images(test_border_masks[0:20,:,:,:],5),'borders')#.show()
# visualize(group_images(img_truth[0:20,:,:,:],5),'gtruth')#.show()#============ Load the data and divide in patches
#图像分块
patches_imgs_test = None
new_height = None
new_width = None
masks_test  = None
patches_masks_test = None
if average_mode == True:patches_imgs_test, new_height, new_width, masks_test = get_data_testing_overlap(DRIVE_test_imgs_original=DRIVE_test_imgs_original,  #originalDRIVE_test_groudTruth=path_data + config.get('data paths', 'test_groundTruth'),  #masksImgs_to_test=int(config.get('testing settings', 'full_images_to_test')),patch_height=patch_height,patch_width=patch_width,stride_height=stride_height,stride_width=stride_width)
else:patches_imgs_test, patches_masks_test = get_data_testing(DRIVE_test_imgs_original=DRIVE_test_imgs_original,  #originalDRIVE_test_groudTruth=path_data + config.get('data paths', 'test_groundTruth'),  #masksImgs_to_test=int(config.get('testing settings', 'full_images_to_test')),patch_height=patch_height,patch_width=patch_width,)#================ Run the prediction of the patches ==================================
best_last = config.get('testing settings', 'best_last')
#Load the saved model
#加载已经训练好的模型 和 相关的权重
model = model_from_json(open(path_experiment+name_experiment +'_architecture.json').read())
model.load_weights(path_experiment+name_experiment + '_'+best_last+'_weights.h5')
#Calculate the predictions
#进行模型预测
predictions = model.predict(patches_imgs_test, batch_size=32, verbose=2)
# verbose = 1 采用进度条形式进行显示  verbose = 2 为每个epoch输出一行记录
print("predicted images size :")
print(predictions.shape)#===== Convert the prediction arrays in corresponding images
pred_patches = pred_to_imgs(predictions, patch_height, patch_width, "original")#========== Elaborate and visualize the predicted images ====================
pred_imgs = None
orig_imgs = None
gtruth_masks = None
if average_mode == True:pred_imgs = recompone_overlap(pred_patches, new_height, new_width, stride_height, stride_width)# predictionsorig_imgs = my_PreProc(test_imgs_orig[0:pred_imgs.shape[0],:,:,:])    #originalsgtruth_masks = masks_test  #ground truth masks
else:pred_imgs = recompone(pred_patches,13,12)       # predictionsorig_imgs = recompone(patches_imgs_test,13,12)  # originalsgtruth_masks = recompone(patches_masks_test,13,12)  #masks
# apply the DRIVE masks on the repdictions #set everything outside the FOV to zero!!
kill_border(pred_imgs, test_border_masks)  #DRIVE MASK  #only for visualization
## back to original dimensions
orig_imgs = orig_imgs[:,:,0:full_img_height,0:full_img_width]
pred_imgs = pred_imgs[:,:,0:full_img_height,0:full_img_width]
gtruth_masks = gtruth_masks[:,:,0:full_img_height,0:full_img_width]
print("Orig imgs shape: " +str(orig_imgs.shape))
print("pred imgs shape: " +str(pred_imgs.shape))
print("Gtruth imgs shape: " +str(gtruth_masks.shape))
visualize(group_images(orig_imgs,N_visual),path_experiment+"all_originals")#.show()
visualize(group_images(pred_imgs,N_visual),path_experiment+"all_predictions")#.show()
visualize(group_images(gtruth_masks,N_visual),path_experiment+"all_groundTruths")#.show()
#visualize results comparing mask and prediction:
#可视化结果 对比预测 与 金标准
assert (orig_imgs.shape[0]==pred_imgs.shape[0] and orig_imgs.shape[0]==gtruth_masks.shape[0])
N_predicted = orig_imgs.shape[0]
group = N_visual
assert (N_predicted%group==0)
for i in range(int(N_predicted/group)):orig_stripe = group_images(orig_imgs[i*group:(i*group)+group,:,:,:],group)masks_stripe = group_images(gtruth_masks[i*group:(i*group)+group,:,:,:],group)pred_stripe = group_images(pred_imgs[i*group:(i*group)+group,:,:,:],group)total_img = np.concatenate((orig_stripe,masks_stripe,pred_stripe),axis=0)visualize(total_img,path_experiment+name_experiment +"_Original_GroundTruth_Prediction"+str(i))#.show()#====== Evaluate the results
#主要用了sklearn模块的中模型评价函数， sklearn.metrics
print("\n\n=======================  Evaluate the results =======================")
#predictions only inside the FOV
#只预测FOV内的部分
y_scores, y_true = pred_only_FOV(pred_imgs,gtruth_masks, test_border_masks)  #returns data only inside the FOV
print("Calculating results only inside the FOV:")
print("y scores pixels: " +str(y_scores.shape[0]) +" (radius 270: 270*270*3.14==228906), including background around retina: " +str(pred_imgs.shape[0]*pred_imgs.shape[2]*pred_imgs.shape[3]) +" (584*565==329960)")
print("y true pixels: " +str(y_true.shape[0]) +" (radius 270: 270*270*3.14==228906), including background around retina: " +str(gtruth_masks.shape[2]*gtruth_masks.shape[3]*gtruth_masks.shape[0])+" (584*565==329960)")#Area under the ROC curve
fpr, tpr, thresholds = roc_curve((y_true), y_scores)
AUC_ROC = roc_auc_score(y_true, y_scores)
# test_integral = np.trapz(tpr,fpr) #trapz is numpy integration
print("\nArea under the ROC curve: " +str(AUC_ROC))
roc_curve =plt.figure()
plt.plot(fpr,tpr,'-',label='Area Under the Curve (AUC = %0.4f)' % AUC_ROC)
plt.title('ROC curve')
plt.xlabel("FPR (False Positive Rate)")
plt.ylabel("TPR (True Positive Rate)")
plt.legend(loc="lower right")
plt.savefig(path_experiment+"ROC.png")#Precision-recall curve
#sklearn.metrics.precision_recall_curve：精确度-召回率曲线
precision, recall, thresholds = precision_recall_curve(y_true, y_scores)
precision = np.fliplr([precision])[0]  #so the array is increasing (you won't get negative AUC)
recall = np.fliplr([recall])[0]  #so the array is increasing (you won't get negative AUC)
AUC_prec_rec = np.trapz(precision,recall)
print("\nArea under Precision-Recall curve: " +str(AUC_prec_rec))
prec_rec_curve = plt.figure()
plt.plot(recall,precision,'-',label='Area Under the Curve (AUC = %0.4f)' % AUC_prec_rec)
plt.title('Precision - Recall curve')
plt.xlabel("Recall")
plt.ylabel("Precision")
plt.legend(loc="lower right")
plt.savefig(path_experiment+"Precision_recall.png")#Confusion matrix
#sklearn.metrics.confusion_matrix ： 混淆矩阵
#混淆矩阵是对有监督学习分类算法准确率进行评估的工具。通过将模型预测的数据与测试数据进行对比，使用各种指标对模型的分类效果进行度量。
threshold_confusion = 0.5
print("\nConfusion matrix:  Costum threshold (for positive) of " +str(threshold_confusion))
y_pred = np.empty((y_scores.shape[0]))
for i in range(y_scores.shape[0]):if y_scores[i]>=threshold_confusion:y_pred[i]=1else:y_pred[i]=0
confusion = confusion_matrix(y_true, y_pred)
print(confusion)
accuracy = 0
if float(np.sum(confusion))!=0:accuracy = float(confusion[0,0]+confusion[1,1])/float(np.sum(confusion))
print("Global Accuracy: " +str(accuracy))
specificity = 0
if float(confusion[0,0]+confusion[0,1])!=0:specificity = float(confusion[0,0])/float(confusion[0,0]+confusion[0,1])
print("Specificity: " +str(specificity))
sensitivity = 0
if float(confusion[1,1]+confusion[1,0])!=0:sensitivity = float(confusion[1,1])/float(confusion[1,1]+confusion[1,0])
print("Sensitivity: " +str(sensitivity))
precision = 0
if float(confusion[1,1]+confusion[0,1])!=0:precision = float(confusion[1,1])/float(confusion[1,1]+confusion[0,1])
print("Precision: " +str(precision))#Jaccard similarity index
#sklearn.metrics. jaccard_similarity_score : jacaard相似度
#jaccard index又称为jaccard similarity coefficient用于比较有限样本集之间的相似性和差异性。给定两个集合A,B jaccard 系数定义为A与B交集的大小与并集大小的比值。jaccard值越大说明相似度越高。jaccard_index = jaccard_similarity_score(y_true, y_pred, normalize=True)
print("\nJaccard similarity score: " +str(jaccard_index))#F1 score F1-score:  是准确率与召回率的综合。 可以认为是平均效果
F1_score = f1_score(y_true, y_pred, labels=None, average='binary', sample_weight=None)
print("\nF1 score (F-measure): " +str(F1_score))#Save the results 保存数据结果
file_perf = open(path_experiment+'performances.txt', 'w')
file_perf.write("Area under the ROC curve: "+str(AUC_ROC)+ "\nArea under Precision-Recall curve: " +str(AUC_prec_rec)+ "\nJaccard similarity score: " +str(jaccard_index)+ "\nF1 score (F-measure): " +str(F1_score)+"\n\nConfusion matrix:"+str(confusion)+"\nACCURACY: " +str(accuracy)+"\nSENSITIVITY: " +str(sensitivity)+"\nSPECIFICITY: " +str(specificity)+"\nPRECISION: " +str(precision))
file_perf.close()

help_functions.py

import h5py
import numpy as np
from PIL import Image
from matplotlib import pyplot as pltdef load_hdf5(infile):        #加载hdf5文件with h5py.File(infile,"r") as f:  #"with" close the file after its nested commands # "image"是写入的时候规定的字段 key-valuereturn f["image"][()]           # 调用方法 train_imgs_original = load_hdf5( file_dir )def write_hdf5(arr,outfile):  #写入hdf5文件with h5py.File(outfile,"w") as f:f.create_dataset("image", data=arr, dtype=arr.dtype)#convert RGB image in black and white
# 将RGB图像转换为Gray图像
def rgb2gray(rgb):assert (len(rgb.shape)==4)  #4D arrays #[Nimgs, channels, height, width]assert (rgb.shape[1]==3)    #确定是否为RGB图像bn_imgs = rgb[:,0,:,:]*0.299 + rgb[:,1,:,:]*0.587 + rgb[:,2,:,:]*0.114bn_imgs = np.reshape(bn_imgs,(rgb.shape[0],1,rgb.shape[2],rgb.shape[3])) # 确保张量形式正确return bn_imgs#group a set of images row per columns
#利用已知信息进行分组显示
#对数据集划分，进行分组显示，totimg图像阵列
def group_images(data,per_row):  # data：数据  per_row：每行显示的图像个数assert data.shape[0]%per_row==0 # data=[Nimgs, channels, height, width]assert (data.shape[1]==1 or data.shape[1]==3)data = np.transpose(data,(0,2,3,1))  #corect format for imshow   # 用于显示all_stripe = []for i in range(int(data.shape[0]/per_row)): # data.shape[0]/per_row 行数stripe = data[i*per_row]  # 相当于matlab中的 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)  # 每per_row张图像拼成一行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) # 每行图像进行拼凑 共len(all_stripe)行return totimg#visualize image (as PIL image, NOT as matplotlib!)
def visualize(data,filename):assert (len(data.shape)==3) #height*width*channelsimg = Noneif data.shape[2]==1:  #in case it is black and whitedata = 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-255else:img = Image.fromarray((data*255).astype(np.uint8))  #the image is between 0-1img.save(filename + '.png') #保存return img#prepare the mask in the right shape for the Unet
# 将金标准图像改写成模型输出形式
def masks_Unet(masks): # size=[Npatches, 1, patch_height, patch_width]assert (len(masks.shape)==4)  #4D arraysassert (masks.shape[1]==1 )  #check the channel is 1im_h = masks.shape[2]im_w = masks.shape[3]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]=0new_masks[i,j,1]=1return new_masks# 网络输出转换成图像子块
# 网络输出 size=[Npatches, patch_height*patch_width, 2]
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] #pred[:, :, 0] 是反分割图像输出 pred[:, :, 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]=1else:pred_images[i,pix]=0else:print("mode " +str(mode) +" not recognized, it can be 'original' or 'threshold'")exit()# 改写成(Npatches,1, height, width)pred_images = np.reshape(pred_images,(pred_images.shape[0],1, patch_height, patch_width))return pred_images

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[1]==3)  #Use the original images#black-white conversiontrain_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 rangereturn train_imgs#============================================================
#========= PRE PROCESSING FUNCTIONS ========================#
#============================================================#==== histogram equalization
def histo_equalized(imgs):assert (len(imgs.shape)==4)  #4D arraysassert (imgs.shape[1]==1)  #check the channel is 1imgs_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 arraysassert (imgs.shape[1]==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 arraysassert (imgs.shape[1]==1)  #check the channel is 1imgs_normalized = np.empty(imgs.shape)imgs_std = np.std(imgs)imgs_mean = np.mean(imgs)imgs_normalized = (imgs-imgs_mean)/imgs_stdfor 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])))*255return imgs_normalizeddef adjust_gamma(imgs, gamma=1.0):assert (len(imgs.shape)==4)  #4D arraysassert (imgs.shape[1]==1)  #check the channel is 1# build a lookup table mapping the pixel values [0, 255] to# their adjusted gamma valuesinvGamma = 1.0 / gammatable = np.array([((i / 255.0) ** invGamma) * 255 for i in np.arange(0, 256)]).astype("uint8")# apply gamma correction using the lookup tablenew_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

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 traintrain_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  # 图像裁剪 size=565*565train_masks = train_masks[:,:,9:574,:]  #cut bottom and top so now it is 565*565 # 图像裁剪 size=565*565data_consistency_check(train_imgs,train_masks) # 训练图像和金标准图像一致性检查#check masks are within 0-1assert(np.min(train_masks)==0 and np.max(train_masks)==1) #金标准图像 2类 0-1print("\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#Load the original data and return the extracted patches for training/testing
def get_data_testing(DRIVE_test_imgs_original, DRIVE_test_groudTruth, Imgs_to_test, patch_height, patch_width):### testtest_imgs_original = load_hdf5(DRIVE_test_imgs_original)test_masks = load_hdf5(DRIVE_test_groudTruth)test_imgs = my_PreProc(test_imgs_original)test_masks = test_masks/255.#extend both images and masks so they can be divided exactly by the patches dimensionstest_imgs = test_imgs[0:Imgs_to_test,:,:,:]test_masks = test_masks[0:Imgs_to_test,:,:,:]test_imgs = paint_border(test_imgs,patch_height,patch_width)test_masks = paint_border(test_masks,patch_height,patch_width)data_consistency_check(test_imgs, test_masks)#check masks are within 0-1assert(np.max(test_masks)==1  and np.min(test_masks)==0)print("\ntest images/masks shape:")print(test_imgs.shape)print("test images range (min-max): " +str(np.min(test_imgs)) +' - '+str(np.max(test_imgs)))print("test masks are within 0-1\n")#extract the TEST patches from the full imagespatches_imgs_test = extract_ordered(test_imgs,patch_height,patch_width)patches_masks_test = extract_ordered(test_masks,patch_height,patch_width)data_consistency_check(patches_imgs_test, patches_masks_test)print("\ntest PATCHES images/masks shape:")print(patches_imgs_test.shape)print("test PATCHES images range (min-max): " +str(np.min(patches_imgs_test)) +' - '+str(np.max(patches_imgs_test)))return patches_imgs_test, patches_masks_test# Load the original data and return the extracted patches for testing
# return the ground truth in its original shape
#测试图像分块
def get_data_testing_overlap(DRIVE_test_imgs_original, DRIVE_test_groudTruth, Imgs_to_test, patch_height, patch_width, stride_height, stride_width):### testtest_imgs_original = load_hdf5(DRIVE_test_imgs_original)test_masks = load_hdf5(DRIVE_test_groudTruth)test_imgs = my_PreProc(test_imgs_original)test_masks = test_masks/255.#extend both images and masks so they can be divided exactly by the patches dimensionstest_imgs = test_imgs[0:Imgs_to_test,:,:,:]test_masks = test_masks[0:Imgs_to_test,:,:,:]test_imgs = paint_border_overlap(test_imgs, patch_height, patch_width, stride_height, stride_width) # 拓展图像 可以准确划分#check masks are within 0-1assert(np.max(test_masks)==1  and np.min(test_masks)==0)print("\ntest images shape:")print(test_imgs.shape)print("\ntest mask shape:")print(test_masks.shape)print("test images range (min-max): " +str(np.min(test_imgs)) +' - '+str(np.max(test_imgs)))print("test masks are within 0-1\n")#extract the TEST patches from the full images# 按照顺序提取图像快 方便后续进行图像恢复（作者采用了overlap策略）patches_imgs_test = extract_ordered_overlap(test_imgs,patch_height,patch_width,stride_height,stride_width)print("\ntest PATCHES images shape:")print(patches_imgs_test.shape)print("test PATCHES images range (min-max): " +str(np.min(patches_imgs_test)) +' - '+str(np.max(patches_imgs_test)))return patches_imgs_test, test_imgs.shape[2], test_imgs.shape[3], test_masks #原始大小#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[2]==masks.shape[2])assert(imgs.shape[3]==masks.shape[3])assert(masks.shape[1]==1)assert(imgs.shape[1]==1 or imgs.shape[1]==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[1]==1 or full_imgs.shape[1]==3)  #check the channel is 1 or 3  # 通道检验assert (full_masks.shape[1]==1)   #masks only black and white  # 通道检验assert (full_imgs.shape[2] == full_masks.shape[2] and full_imgs.shape[3] == full_masks.shape[3])  # 尺寸检验patches = np.empty((N_patches,full_imgs.shape[1],patch_h,patch_w))  # 训练图像总子块patches_masks = np.empty((N_patches,full_masks.shape[1],patch_h,patch_w))  # 训练金标准总子块img_h = full_imgs.shape[2]  #height of the full imageimg_w = full_imgs.shape[3] #width of the full image# (0,0) in the center of the imagepatch_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 FOVif inside==True:if is_patch_inside_FOV(x_center,y_center,img_w,img_h,patch_h)==False:continuepatch = 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  # size=[Npatches, 3, patch_h, patch_w]patches_masks[iter_tot]=patch_mask  # size=[Npatches, 1, patch_h, patch_w]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 centery_ = y - int(img_h/2)  # origin (0,0) shifted to image centerR_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 FOVradius = np.sqrt((x_*x_)+(y_*y_))if radius < R_inside:return Trueelse:return False#Divide all the full_imgs in pacthes
def extract_ordered(full_imgs, patch_h, patch_w):assert (len(full_imgs.shape)==4)  #4D arraysassert (full_imgs.shape[1]==1 or full_imgs.shape[1]==3)  #check the channel is 1 or 3img_h = full_imgs.shape[2]  #height of the full imageimg_w = full_imgs.shape[3] #width of the full imageN_patches_h = int(img_h/patch_h) #round to lowest intif (img_h%patch_h != 0):print("warning: " +str(N_patches_h) +" patches in height, with about " +str(img_h%patch_h) +" pixels left over")N_patches_w = int(img_w/patch_w) #round to lowest intif (img_h%patch_h != 0):print("warning: " +str(N_patches_w) +" patches in width, with about " +str(img_w%patch_w) +" pixels left over")print("number of patches per image: " +str(N_patches_h*N_patches_w))N_patches_tot = (N_patches_h*N_patches_w)*full_imgs.shape[0]patches = np.empty((N_patches_tot,full_imgs.shape[1],patch_h,patch_w))iter_tot = 0   #iter over the total number of patches (N_patches)for i in range(full_imgs.shape[0]):  #loop over the full imagesfor h in range(N_patches_h):for w in range(N_patches_w):patch = full_imgs[i,:,h*patch_h:(h*patch_h)+patch_h,w*patch_w:(w*patch_w)+patch_w]patches[iter_tot]=patchiter_tot +=1   #totalassert (iter_tot==N_patches_tot)return patches  #array with all the full_imgs divided in patches#原始图像进行拓展填充：
def paint_border_overlap(full_imgs, patch_h, patch_w, stride_h, stride_w):assert (len(full_imgs.shape)==4)  #4D arraysassert (full_imgs.shape[1]==1 or full_imgs.shape[1]==3)  #check the channel is 1 or 3img_h = full_imgs.shape[2]  #height of the full imageimg_w = full_imgs.shape[3] #width of the full imageleftover_h = (img_h-patch_h)%stride_h  #leftover on the h dimleftover_w = (img_w-patch_w)%stride_w  #leftover on the w dimif (leftover_h != 0):  #change dimension of img_hprint("\nthe side H is not compatible with the selected stride of " +str(stride_h))print("img_h " +str(img_h) + ", patch_h " +str(patch_h) + ", stride_h " +str(stride_h))print("(img_h - patch_h) MOD stride_h: " +str(leftover_h))print("So the H dim will be padded with additional " +str(stride_h - leftover_h) + " pixels")tmp_full_imgs = np.zeros((full_imgs.shape[0],full_imgs.shape[1],img_h+(stride_h-leftover_h),img_w))tmp_full_imgs[0:full_imgs.shape[0],0:full_imgs.shape[1],0:img_h,0:img_w] = full_imgsfull_imgs = tmp_full_imgsif (leftover_w != 0):   #change dimension of img_wprint("the side W is not compatible with the selected stride of " +str(stride_w))print("img_w " +str(img_w) + ", patch_w " +str(patch_w) + ", stride_w " +str(stride_w))print("(img_w - patch_w) MOD stride_w: " +str(leftover_w))print("So the W dim will be padded with additional " +str(stride_w - leftover_w) + " pixels")tmp_full_imgs = np.zeros((full_imgs.shape[0],full_imgs.shape[1],full_imgs.shape[2],img_w+(stride_w - leftover_w)))tmp_full_imgs[0:full_imgs.shape[0],0:full_imgs.shape[1],0:full_imgs.shape[2],0:img_w] = full_imgsfull_imgs = tmp_full_imgsprint("new full images shape: \n" +str(full_imgs.shape))return full_imgs#Divide all the full_imgs in pacthes
#按顺序提取图像子块：
# 按照顺序对拓展后的图像进行子块采样
def extract_ordered_overlap(full_imgs, patch_h, patch_w,stride_h,stride_w):assert (len(full_imgs.shape)==4)  #4D arraysassert (full_imgs.shape[1]==1 or full_imgs.shape[1]==3)  #check the channel is 1 or 3img_h = full_imgs.shape[2]  #height of the full imageimg_w = full_imgs.shape[3] #width of the full imageassert ((img_h-patch_h)%stride_h==0 and (img_w-patch_w)%stride_w==0)N_patches_img = ((img_h-patch_h)//stride_h+1)*((img_w-patch_w)//stride_w+1)  #// --> division between integers  # 每张图像采集到的子图像N_patches_tot = N_patches_img*full_imgs.shape[0]  # 测试集总共的子图像数量print("Number of patches on h : " +str(((img_h-patch_h)//stride_h+1)))print("Number of patches on w : " +str(((img_w-patch_w)//stride_w+1)))print("number of patches per image: " +str(N_patches_img) +", totally for this dataset: " +str(N_patches_tot))patches = np.empty((N_patches_tot,full_imgs.shape[1],patch_h,patch_w))iter_tot = 0   #iter over the total number of patches (N_patches)for i in range(full_imgs.shape[0]):  #loop over the full imagesfor h in range((img_h-patch_h)//stride_h+1):for w in range((img_w-patch_w)//stride_w+1):patch = full_imgs[i,:,h*stride_h:(h*stride_h)+patch_h,w*stride_w:(w*stride_w)+patch_w]patches[iter_tot]=patchiter_tot +=1   #totalassert (iter_tot==N_patches_tot)return patches  #array with all the full_imgs divided in patches#对于图像子块进行复原
# [Npatches, 1, patch_h, patch_w]  img_h=new_height[588] img_w=new_width[568] stride-[10,10]
def recompone_overlap(preds, img_h, img_w, stride_h, stride_w):assert (len(preds.shape)==4)  #4D arrays # 检查张量尺寸assert (preds.shape[1]==1 or preds.shape[1]==3)  #check the channel is 1 or 3patch_h = preds.shape[2]patch_w = preds.shape[3]N_patches_h = (img_h-patch_h)//stride_h+1 # img_h方向包括的patch_h数量N_patches_w = (img_w-patch_w)//stride_w+1 # img_w方向包括的patch_w数量N_patches_img = N_patches_h * N_patches_w # 每张图像包含的patch的数目print("N_patches_h: " +str(N_patches_h))print("N_patches_w: " +str(N_patches_w))print("N_patches_img: " +str(N_patches_img))assert (preds.shape[0]%N_patches_img==0)N_full_imgs = preds.shape[0]//N_patches_img  # 全幅图像的数目print("According to the dimension inserted, there are " +str(N_full_imgs) +" full images (of " +str(img_h)+"x" +str(img_w) +" each)")full_prob = np.zeros((N_full_imgs,preds.shape[1],img_h,img_w))  #itialize to zero mega array with sum of Probabilitiesfull_sum = np.zeros((N_full_imgs,preds.shape[1],img_h,img_w))k = 0 #iterator over all the patches  #迭代所有的子块for i in range(N_full_imgs):for h in range((img_h-patch_h)//stride_h+1):for w in range((img_w-patch_w)//stride_w+1):full_prob[i,:,h*stride_h:(h*stride_h)+patch_h,w*stride_w:(w*stride_w)+patch_w]+=preds[k]full_sum[i,:,h*stride_h:(h*stride_h)+patch_h,w*stride_w:(w*stride_w)+patch_w]+=1k+=1assert(k==preds.shape[0])assert(np.min(full_sum)>=1.0)  #at least onefinal_avg = full_prob/full_sum # 叠加概率 / 叠加权重 ： 采用了均值的方法print(final_avg.shape)assert(np.max(final_avg)<=1.0) #max value for a pixel is 1.0assert(np.min(final_avg)>=0.0) #min value for a pixel is 0.0return final_avg#Recompone the full images with the patches
def recompone(data,N_h,N_w):assert (data.shape[1]==1 or data.shape[1]==3)  #check the channel is 1 or 3assert(len(data.shape)==4)N_pacth_per_img = N_w*N_hassert(data.shape[0]%N_pacth_per_img == 0)N_full_imgs = data.shape[0]/N_pacth_per_imgpatch_h = data.shape[2]patch_w = data.shape[3]N_pacth_per_img = N_w*N_h#define and start full recomponefull_recomp = np.empty((N_full_imgs,data.shape[1],N_h*patch_h,N_w*patch_w))k = 0  #iter full imgs = 0  #iter single patchwhile (s<data.shape[0]):#recompone one:single_recon = np.empty((data.shape[1],N_h*patch_h,N_w*patch_w))for h in range(N_h):for w in range(N_w):single_recon[:,h*patch_h:(h*patch_h)+patch_h,w*patch_w:(w*patch_w)+patch_w]=data[s]s+=1full_recomp[k]=single_reconk+=1assert (k==N_full_imgs)return full_recomp#Extend the full images because patch divison is not exact
def paint_border(data,patch_h,patch_w):assert (len(data.shape)==4)  #4D arraysassert (data.shape[1]==1 or data.shape[1]==3)  #check the channel is 1 or 3img_h=data.shape[2]img_w=data.shape[3]new_img_h = 0new_img_w = 0if (img_h%patch_h)==0:new_img_h = img_helse:new_img_h = ((int(img_h)/int(patch_h))+1)*patch_hif (img_w%patch_w)==0:new_img_w = img_welse:new_img_w = ((int(img_w)/int(patch_w))+1)*patch_wnew_data = np.zeros((data.shape[0],data.shape[1],new_img_h,new_img_w))new_data[:,:,0:img_h,0:img_w] = data[:,:,:,:]return new_data#return only the pixels contained in the FOV, for both images and masks
def pred_only_FOV(data_imgs,data_masks,original_imgs_border_masks):assert (len(data_imgs.shape)==4 and len(data_masks.shape)==4)  #4D arraysassert (data_imgs.shape[0]==data_masks.shape[0])assert (data_imgs.shape[2]==data_masks.shape[2])assert (data_imgs.shape[3]==data_masks.shape[3])assert (data_imgs.shape[1]==1 and data_masks.shape[1]==1)  #check the channel is 1height = data_imgs.shape[2]width = data_imgs.shape[3]new_pred_imgs = []new_pred_masks = []for i in range(data_imgs.shape[0]):  #loop over the full imagesfor x in range(width):for y in range(height):if inside_FOV_DRIVE(i,x,y,original_imgs_border_masks)==True:new_pred_imgs.append(data_imgs[i,:,y,x])new_pred_masks.append(data_masks[i,:,y,x])new_pred_imgs = np.asarray(new_pred_imgs)new_pred_masks = np.asarray(new_pred_masks)return new_pred_imgs, new_pred_masks#function to set to black everything outside the FOV, in a full image
def kill_border(data, original_imgs_border_masks):assert (len(data.shape)==4)  #4D arraysassert (data.shape[1]==1 or data.shape[1]==3)  #check the channel is 1 or 3height = data.shape[2]width = data.shape[3]for i in range(data.shape[0]):  #loop over the full imagesfor x in range(width):for y in range(height):if inside_FOV_DRIVE(i,x,y,original_imgs_border_masks)==False:data[i,:,y,x]=0.0def inside_FOV_DRIVE(i, x, y, DRIVE_masks):assert (len(DRIVE_masks.shape)==4)  #4D arraysassert (DRIVE_masks.shape[1]==1)  #DRIVE masks is black and white# DRIVE_masks = DRIVE_masks/255.  #NOOO!! otherwise with float numbers takes forever!!if (x >= DRIVE_masks.shape[3] or y >= DRIVE_masks.shape[2]): #my image bigger than the originalreturn Falseif (DRIVE_masks[i,0,y,x]>0):  #0==black pixels# print DRIVE_masks[i,0,y,x]  #verify it is working rightreturn Trueelse:return False

Unet实现图像分割（四）相关推荐

使用U-Net 进行图像分割
最近做病理AI的细胞计数问题,需要对图像中的各个细胞进行分类,若采用普通的CNN+普通图像分割,估计实现效果不佳.为了解决这个问题,大致有两种方案:目标检测和图像分割.目标检测的算法以Faster ...
keras搭建Unet实现图像分割
本文运行系统为Windows,若使用Linux请注意相关指令及操作. 目录一.环境配置 1.创建及激活虚拟环境. 2.安装相关库二.Unet原理讲解 1.结构简析 2.结构细述三.模型训练及测试 ...
nnU-Net: 基于U-Net医学图像分割技术的自适应框架
** nnU-Net: 基于U-Net医学图像分割技术的自适应框架 ** https://arxiv.org/pdf/1809.10486.pdf 作者:Fabian Isensee 提要 U-Net ...
TensorFlow实现Unet遥感图像分割
Unet是一种U型网络,分为左右两部分卷积,左边为下采样提取高维特征,右边为上采样并与左侧融合实现图像分割.这里使用TensorFlow实现Unet网络,实现对遥感影像的道路分割. 训练数据: 标签图 ...
Tensorflow2.0 U-Net医学图像分割(胸部X光图像分割)
1.U-Net网络结构论文题目:<U-Net: Convolutional Networks for Biomedical Image Segmentation> 论文地址:https: ...
Unet实现图像分割（一）
在Rrtina-Unet-master文件夹下新建一个test文件夹(源代码没有这个test文件夹会报错) 将lib文件夹下的三个py文件里面的from xxx import xxx格式改成一致,就是 ...
U-net进行图像分割
数据准备原始数据:首先准备数据,参考数据来自于 ISBI 挑战的数据集.数据可以在这里下载到,含30张训练图.30张对应的标签.30张测试图片,均为.tif 格式文件程序准备程序地址:这里 ...
Unet实现图像分割（三）
1.retinaNN_training.py 模型Model的compile方法: model.compile(self, optimizer, loss, metrics=None, loss_we ...
Unet实现图像分割（二）
使用了google colaboratory的免费GPU进行训练,调整了源代码的各个参数,下面是configuration.txt文件的解析: [data paths] 只有在修改了prepare_d ...

Unet实现图像分割（四）

Unet实现图像分割（四）相关推荐

最新文章

热门文章