Kaggle肾小球图像分割比赛全解析

概览

赛事描述

Kaggle上正在进行一项名为“HuBMAP: Hacking the Kidney”的竞赛,链接见这里,总奖金是6万美金,该竞赛的目的是为了检测人体组织中的functional tissue units (FTUs)。FTU的定义是:three-dimensional block of cells centered around a capillary, such that each cell in this block is within diffusion distance from any other cell in the same block,感觉类似于细胞的概念,正中是细胞核,周围是细胞质。

算法评估标准

该竞赛使用Dice系数来评估算法的优劣。关于Dice系数,可以见如下博客解析:
医学图像分割之 Dice Loss

竞赛所提交的文件使用游程编码方式(RLE,run-length encoding)来减小文件体积。
关于掩膜mask与rle编码与解码的代码,可以参见网友Paulo Pinto的notebook:
RLE functions - Run Lenght Encode & Decode

数据分析

数据集概览

该赛事中的数据一共有20张肾的图像,每一张都对其中的肾小球FTU进行了标注,有8张用于训练集,5张用于公榜测试集,剩下7张用于私榜测试集。
每一张图像都是非常大的TIFF格式,500MB-5GB大小。

训练集中的标注有两种形式:游程编码和未编码的JSON格式。
可以使用外部数据和/或预训练的机器学习模型,不过这些数据和模型必须在CC BY 4.0下授权。
下载数据集(这是在google colab上运行,colab上的机器性能较好;如果直接使用kaggle上的notebook,则数据集直接内置)

1
!kaggle competitions download -c hubmap-kidney-segmentation

然后进行数据的探索性分析,该过程参考了以下notebook:
HuBMAP - Exploratory Data Analysis

导入必要的包

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
import pandas as pd
import pathlib, sys, os, random, time
import numba, cv2, gc

import matplotlib.pyplot as plt
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

from tqdm.notebook import tqdm

import albumentations as A
import rasterio
from rasterio.windows import Window

import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.utils.data as D

import torchvision
from torchvision import transforms as T

配置路径及超参数

1
2
3
4
5
6
7
BASE_PATH = "../input/hubmap-kidney-segmentation/"
TRAIN_PATH = os.path.join(BASE_PATH, "train/")
TEST_PATH = os.path.join(BASE_PATH, "test/")

EPOCHES = 5
BATCH_SIZE = 32
DEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'

文件分析

(1)训练集文件分析:

1
2
df_train = pd.read_csv(os.path.join(BASE_PATH, "train.csv"))
df_train

得到:

1
2
3
4
5
6
7
8
9
           id                        encoding
0 2f6ecfcdf 296084587 4 296115835 6 296115859 14 296147109...
1 aaa6a05cc 30989109 59 31007591 64 31026074 68 31044556 7...
2 cb2d976f4 78144363 5 78179297 15 78214231 25 78249165 35...
3 0486052bb 101676003 6 101701785 8 101727568 9 101753351 ...
4 e79de561c 7464094 14 7480273 41 7496453 67 7512632 82 75...
5 095bf7a1f 113430380 22 113468538 67 113506697 111 113544...
6 54f2eec69 124601765 36 124632133 109 124662536 147 12469...
7 1e2425f28 49453112 7 49479881 22 49506657 31 49533433 40...

train.csv文件中包含了图像的id及其游程编码。可以看出图像的名称就是id名。
(2)提交文件分析:

1
2
df_sub = pd.read_csv(os.path.join(BASE_PATH, "sample_submission.csv"))
df_sub

得到:

1
2
3
4
5
6
   id                          predicted
0 b9a3865fc NaN
1 b2dc8411c NaN
2 26dc41664 NaN
3 c68fe75ea NaN
4 afa5e8098 NaN

提交文件就是对公共测试集上的图像的预测,可以看出id就是公共测试集中的图像名称,predicted一栏需要后面填入。
(3)数据集大小分析:

1
2
print("number of train images: ", df_train.shape[0])
print("number of test images: ", df_sub.shape[0])

分别是8和5。
(4)元数据分析:

1
2
df_meta = pd.read_csv(os.path.join(BASE_PATH, "HuBMAP-20-dataset_information.csv"))
df_meta.sample(3)

该文件中包含了数据集中的每一张图像额外的信息,比如它的主人的身体信息、性别、种族等。
同时指明训练集中除了肾小球的标注文件,比如1e2425f28.json,还有其他解剖组织的标注文件,比如1e2425f28-anatomical-structure.json。
该文件是为了辅助理解该赛题背后的医学知识,有可能对特征功能有用,但目前看没法直接使用。

工具函数

以下是关于游程编码和解码、读取图像、可视化的工具代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
# 图像分块
# min_overlap这个参数指的是有可能出现的最小的overlap,而不是保证这个overlap一定会出现
# 即自适应产生的overlap肯定会大于该min_overlap
def make_grid(shape, window=256, min_overlap=32):
"""
Return Array of size (N,4), where N - number of tiles,
2nd axis represente slices: x1,x2,y1,y2
"""
x, y = shape
nx = x // (window - min_overlap) + 1
x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)
x1[-1] = x - window
x2 = (x1 + window).clip(0, x)
ny = y // (window - min_overlap) + 1
y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)
y1[-1] = y - window
y2 = (y1 + window).clip(0, y)
slices = np.zeros((nx,ny, 4), dtype=np.int64)

for i in range(nx):
for j in range(ny):
slices[i,j] = x1[i], x2[i], y1[j], y2[j]
return slices.reshape(nx*ny,4)

# 固定随机数,以保证可复现性
def set_seeds(seed=42):
random.seed(seed)
os.environ['PYTHONHASHSEED'] = str(seed)
np.random.seed(seed)
torch.manual_seed(seed)
torch.cuda.manual_seed(seed)
torch.cuda.manual_seed_all(seed)
torch.backends.cudnn.deterministic = True

# 游程编码转为图像掩膜
def rle2mask(mask_rle, shape):
# shape的形状是(width, height), width是通常理解的图像宽度
# 原始的rle编码是str类型,里面的元素成对出现,即start起始像素及length长度
s = mask_rle.split()

# 将s中的start和length分别提取出来,因为它们在原始列表中是成对出现,所以这里对这两种数据都是每隔两个元素提取一次
# 然后再通过python的列表生成式语法另存成numpy数组
# https://www.liaoxuefeng.com/wiki/1016959663602400/1017317609699776
starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]

# 原始csv文件中的游程编码是从绝对位置开始,因此转化为numpy数组时需要减1
starts -= 1
ends = starts + lengths

# 根据原始图像大小建立一个空白图像
img = np.zeros(shape[0]*shape[1], dtype=np.uint8)

# 根据起始像素和结束像素,在该空白图像上创建掩膜
# 1为mask,0为背景
for lo, hi in zip(starts, ends):
img[lo: hi] = 1

# 将img按输入shape变换形状,这里就体现了shape中元素顺序的关键
# 因为RLE编码是先从上到下,然后再从左到右进行的,所以分割时是先满足高度要求,即先按一列一列地来分组
# 因为输入的shape是宽度在前,高度在后,正好reshape就按这个shape来变换形状
# 比如原图如果宽为5,高为2,那么就reshape((5, 2)),即分成5组,每组2个元素
# 然后因为图像存成numpy数组时是行数乘以列数,即转置一下即可
return img.reshape(shape).T

# 图像掩膜转为游程编码
# https://blog.csdn.net/qq_35985044/article/details/104332577
def mask2rle(img):
# 将掩膜按从上到下、从左到右的顺序压平
pixels = img.T.flatten()

# 前后各加一个0作为缓冲区
pixels = np.concatenate([[0], pixels, [0]])

# 以下记录的是掩膜值开始发生变化的位置,使用的方法是将数组错移一位,并与原数组比较
# 这样每一段重复的序列在前后位置都有一个变化的位置记录,前后位置是成对出现的
# +1是为了做位置调整
runs = np.where(pixels[1:] != pixels[:-1])[0] + 1

# 这一步比较抽象,首先记住runs记录了像素值发生变化的位置
# runs[1::2]是从第二个位置开始,每隔两个元素提取一次,即该表达式提取的是“成对出现的前后位置”中的“后”这一位置
# runs[::2]则是从第一个位置开始,每隔两个元素提取一次,即该表达式提取的是“成对出现的前后位置”中的“前”这一位置
# 然后两者相减,并在原来“后”这一位置存储差值,即每个重复序列的长度
runs[1::2] -= runs[::2]

# 将上述元素逐个取出,并且使用空格符连接成字符串
return ' '.join(str(x) for x in runs)

# 使用numba加速
# 可以参考如下nb里的讨论部分
# https://www.kaggle.com/leighplt/pytorch-fcn-resnet50/comments
@numba.njit()
def mask2rle_numba_1d(pixels):
size = len(pixels)
points = []
if pixels[0] == 1: points.append(1)
for i in range(1, size):
if pixels[i] != pixels[i-1]:
if len(points) % 2 == 0:
points.append(i+1)
else:
points.append(i+1 - points[-1])

if pixels[-1] == 1: points.append(size-points[-1]+1)
return points

# 该函数必须得与上面的分开,因为最后的join函数不支持numba
def mask2rle_numba(image):
pixels = image.T.flatten()
points = mask2rle_numba_1d(pixels)
return ' '.join(str(x) for x in points)

# 读取训练集图像
# 对于大型tif图像,推荐使用rasterio来读取,读取速度会提升很多倍
def read_image(image_id, scale=None, verbose=1):
# python3.6引入的f-string,用于格式化字符串
# https://blog.csdn.net/sunxb10/article/details/81036693
image = tifffile.imread(os.path.join(TRAIN_PATH, f"{image_id}.tiff"))

# 数据集中有几张图像是(1, 1, 3, X, Y)这样的格式,需要将其转为正确的(X, Y, 3)格式
if len(image.shape) == 5:
image = image.sequeeze().transpose(1, 2, 0)

# 这个地方用到了pandas对象的布尔数组索引
# https://www.pypandas.cn/docs/user_guide/indexing.html#%E7%B4%A2%E5%BC%95%E7%9A%84%E4%B8%8D%E5%90%8C%E9%80%89%E6%8B%A9
# 特别需要注意的是第二个参数的顺序,这里是图像的宽度在前,高度在后
mask = rle2mask(df_train[df_train["id"] == image_id]['encoding'].values[0], (image.shape[1], image.shape[0]))

# 开启详细显示信息模式
if verbose:
print(f"[{image_id}] Image shape: {image.shape}")
print(f"[{image_id}] Mask shape: {mask.shape}")

# 缩放
if scale:
# 注意opencv的resize函数要求的参数是先宽后高,注意顺序
new_size = (image.shape[1] // scale, image.shape[0] // scale)
image = cv2.resize(image, new_size)
mask = cv2.resize(mask, new_size)

if verbose:
print(f"[{image_id}] Image shape: {image.shape}")
print(f"[{image_id}] Mask shape: {mask.shape}")

return image, mask

# 读取测试集图像
def read_test_image(image_id, scale=None, verbose=1):
image = tifffile.imread(os.path.join(TEST_PATH, f"{image_id}.tiff"))

if len(image.shape) == 5:
image = image.sequeeze().tranpose(1, 2, 0)

if verbose:
print(f"[{image_id}] Image shape: {image.shape}")

if scale:
new_size = (image.shape[1] // scale, image.shape[0] // scale)
image = cv2.resize(image, new_size)

if verbose:
print(f"[{image_id}] Image shape: {image.shape}")

return image

# 可视化图像及其掩膜
# 主要是应用了matplotlib库,其用法可参考如下链接
# https://lijin-thu.github.io/06.%20matplotlib/06.01%20pyplot%20tutorial.html
def plot_image_and_mask(image, mask, image_id):
# 产生一幅图,指定其大小
plt.figure(figsize=(16, 10))

# 创建一行三列的子图
# 这里是第一个子图
plt.subplot(1, 3, 1)
plt.imshow(image)
plt.title(f"Image {image_id}", fontsize=18)

# 第二个子图
plt.subplot(1, 3, 2)
plt.imshow(image)
plt.imshow(mask, cmap="hot", alpha=0.5)
plt.title(f"Image {image_id} + mask", fontsize=18)

# 第三个子图
plt.subplot(1, 3, 3)
plt.imshow(mask, cmap="hot")
plt.title(f"Mask", fontsize=18)

plt.show()

# 可视化图像及其掩膜的一部分
def plot_slice_image_and_mask(image, mask, start_h, end_h, start_w, end_w):
plt.figure(figsize=(16, 5))
sub_image = image[start_h : end_h, start_w : end_w, :]
sub_mask = mask[start_h : end_h, start_w : end_w]
plt.subplot(1, 3, 1)
plt.imshow(sub_image)
plt.axis("off")
plt.subplot(1, 3, 2)
plt.imshow(sub_image)
plt.imshow(sub_mask, cmap="hot", alpha=0.5)
plt.axis("off")
plt.subplot(1, 3, 3)
plt.imshow(sub_mask, cmap="hot")
plt.axis("off")
plt.show()

依据这些工具函数,可以很方便地进行数据集的读取和调用。
以训练集中的某张图像为例,观察其部分原图及其标注:

1
2
3
4
5
6
image_id = "0486052bb"
image, mask = read_image(image_id, 2)
plot_image_and_mask(image, mask, image_id)

plot_slice_image_and_mask(image, mask, 5000, 7500, 2500, 5000)
plot_slice_image_and_mask(image, mask, 5250, 5720, 3500, 4000)

结果如下图:
vis

无网络连接安装依赖

因为这个比赛的notebook不允许使用Internet,因此,如果调用了kaggle默认环境所没有的模块和包时,需要事先自己在线下准备好。
具体做法可以参考该教程
(1)首先打开notebook的网络,然后使用pip下载所依赖的包的whl文件:

1
!pip download [package_name]

记住此时download的顺序以及版本号(这点一定要注意),后面要按该顺序的逆序进行安装。
(2)将这些whl文件上传到kaggle的dataset中;
这一步又涉及怎样将这些文件先下载下来,此时可以参考该教程
一种是直接在右侧的output侧栏中点击下载;
一种是通过命令下载:

1
2
%cd /kaggle/working
from IPython.display import FileLink -> FileLink(r'*name of file*')

注意,在kaggle上上传非whl的安装包时,比如上传的时.tar.gz格式的源码包,kaggle会自动将其解压成文件夹。
因此,需要将此种文件后缀名需要更改为.xyz,然后上传。
上传时如果发现其他dataset仓库已经有相同的软件包,可以直接用那里的,也可以选择include duplicates上传自己的。
(3)在该notebook中引用上面的dataset,然后按逆序安装这些whl。
对于非whl的文件,注意将其copy到本地路径后,再修改为.tar.gz后缀名,然后正常pip install即可。

正确提交一次和跑一遍模型

正确提交一次

这里首先通过正确提交一次,保证提交文件是正确的,否则可能花了很多时间搭建算法和训练,最后却卡在提交上,没法及时提交。
因为Kaggle的submission目前看像是一个玄学,大家都在尝试怎样提交成功,比如下面的讨论:
Scoring error on submission - even with sample_submission.csv
Scoring error again
原因就是如第一个帖子中Julia所说,kaggle官方为了确保测试集不会被大家hack或产生leak,将测试集藏得很深。
想要提交成功,要假设测试集就在那里,且文件的ID千万不能硬编码,要做到“自适应”,同时里面的内容一开始可以全部设为0,但也注意有些竞赛对于数值也有要求,比如下面的帖子中,因为最后的score要计算Spearman系数,所以每一个array中至少要有两个不同的值:
Unable to fix submission scoring error
针对于此例,一个很小的提交如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
DATA_PATH = '../input/hubmap-kidney-segmentation'
identity = rasterio.Affine(1, 0, 0, 0, 1, 0)

p = pathlib.Path(DATA_PATH)
subm = {}

for i, filename in enumerate(p.glob('test/*.tiff')):
dataset = rasterio.open(filename.as_posix(), transform = identity)
preds = np.zeros(dataset.shape, dtype=np.uint8)
subm[i] = {'id':filename.stem, 'predicted': mask2rle_numba(preds)}
del preds
gc.collect();

submission = pd.DataFrame.from_dict(subm, orient='index')
submission.to_csv('submission.csv', index=False)

跑一遍模型

这里实际是用随机权重的模型在测试集上跑一遍,保证整个流程是通路的,然后再进行模型的训练,即“以终为始”。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def get_model():
model = torchvision.models.segmentation.fcn_resnet50(False)

# 将原来的用于多类分割的模型最后一层改为目前的两类分割
model.classifier[4] = nn.Conv2d(512, 1, kernel_size=(1, 1), stride=(1, 1))
return model

# 这一步也非常重要,因为不能有网络连接,如果不事先下载好这些权重文件,notebook运行过程中会联网下载,导致提交不成功
!mkdir -p /root/.cache/torch/hub/checkpoints/
!cp ../input/pytorch-pretrained-models/resnet50-19c8e357.pth /root/.cache/torch/hub/checkpoints/
!cp ../input/pretrain-coco-weights-pytorch/fcn_resnet50_coco-1167a1af.pth /root/.cache/torch/hub/checkpoints/

model = get_model()
model.to(DEVICE);

identity = rasterio.Affine(1, 0, 0, 0, 1, 0)

# 这个变换一定要与训练集的变换相同
trfm = T.Compose([
T.ToPILImage(),
T.Resize(NEW_SIZE),
T.ToTensor(),
T.Normalize([0.625, 0.448, 0.688],
[0.131, 0.177, 0.101]),
])

p = pathlib.Path(DATA_PATH)

subm = {}

## 真正使用过程中,在训练模型后,在这里要加载训练好的模型
# model.load_state_dict(torch.load("../input/models/HuBMAP_model_best.pth"))
model.eval()

# 遍历测试集中的文件
for i, filename in enumerate(p.glob('test/*.tiff')):
# 使用rasterio来打开tiff文件,实测速度会很快
# transform参数是用来进行仿射变换,这里不涉及坐标系变换,可以不用
# 常用用法可见:
# https://theonegis.github.io/geos/%E4%BD%BF%E7%94%A8Rasterio%E8%AF%BB%E5%8F%96%E6%A0%85%E6%A0%BC%E6%95%B0%E6%8D%AE/index.html
dataset = rasterio.open(filename.as_posix(), transform = identity)

# 因为原始图像很大,为防止爆内存,将其分块处理
slices = make_grid(dataset.shape, window=WINDOW, min_overlap=MIN_OVERLAP)
# 创建一个存储预测结果的缓存
preds = np.zeros(dataset.shape, dtype=np.uint8)
# 对分块图像进行遍历
for (x1,x2,y1,y2) in slices:
# read方法将rasterio的数据集格式转为numpy.ndarray
image = dataset.read([1,2,3], window=Window.from_slices((x1,x2),(y1,y2)))
# 改变一下维度次序
image = np.moveaxis(image, 0, -1)
# 与训练集采用同样的tranform变换
image = trfm(image)
with torch.no_grad():
# 将image放入DEVICE中,cpu或gpu
image = image.to(DEVICE)[None]
# 模型对图像进行预测
# 这个地方需要注意的是model的返回值
# torchvision模型库中的classification和segmentation、detection的模型返回值不同
# 比如,用于classification的模型,比如AlexNet,其model(X)返回的就是预测值y的tensor
# 而segmentation的模型,比如fcn_resnet50,它的返回值是一个有序字典,其中的key有out和aux
# 所以下面的代码需要使用out这个key来获得预测值tensor
# 可以参考:
# https://pytorch.org/docs/stable/torchvision/models.html
# https://github.com/pytorch/vision/blob/d0063f3d83beac01e85f3027c4de6499a8985469/torchvision/models/segmentation/fcn.py#L9
# https://colab.research.google.com/github/spmallick/learnopencv/blob/master/PyTorch-Segmentation-torchvision/intro-seg.ipynb#scrollTo=ZsIngeXleQ1H
score = model(image)['out'][0][0]

# 这个得分有可能是负的,为了一致性,将其使用sigmoid激活,转为0到1
score_sigmoid = score.sigmoid().cpu().numpy()
score_sigmoid = cv2.resize(score_sigmoid, (WINDOW, WINDOW))

# 以0.5作为阈值,输出mask
preds[x1:x2,y1:y2] = (score_sigmoid > 0.5).astype(np.uint8)

# 将预测值转为RLE编码
subm[i] = {'id':filename.stem, 'predicted': mask2rle_numba(preds)}
# 删除临时变量
del preds
# 回收内存,防止内部爆掉,这一步非常重要
gc.collect();

submission = pd.DataFrame.from_dict(subm, orient='index')
submission.to_csv('submission.csv', index=False)

制作数据集

因为原始图像太大了,单张甚至达到5G大小,极易爆内存,因此需要将其切分成小图像,才能进行可行的训练。

在线制作数据集

这一节之所以称为“在线制作”,是因为数据集的加载和切分都是在程序运行时才开始进行,与之对应的,下面一节采用的是事先将原来的数据集切分好,等到用的时候直接读入即可。
这一节的在线制作是直接制作了适用于PyTorch的数据集格式,如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
identity = rasterio.Affine(1, 0, 0, 0, 1, 0)

# 自定义数据集的制作可以参考如下链接:
# https://pytorch.org/tutorials/beginner/data_loading_tutorial.html

class HubDataset(D.Dataset):
def __init__(self, root_dir, transform,
window=256, overlap=32, threshold = 100):
# 加载路径使用pathlib,后续推荐使用它来替代os.path
# 具体使用方法可以参考:
# https://docs.python.org/zh-cn/3/library/pathlib.html
self.path = pathlib.Path(root_dir)
self.overlap = overlap
self.window = window
self.transform = transform
self.csv = pd.read_csv((self.path / 'train.csv').as_posix(),
index_col=[0])
self.threshold = threshold
self.x, self.y = [], []

# 运行下面的分块函数
self.build_slices()
self.len = len(self.x)
self.as_tensor = T.Compose([
# ToTensor函数接受PIL Image或numpy.ndarray,将其先由HWC转置为CHW格式,再转为float后每个像素除以255,
# 从而将数据范围变换到[0, 1]之间
# https://www.cnblogs.com/ocean1100/p/9494640.html
T.ToTensor(),

# 标准化的操作是减去均值并除以标准差,即将数据变换为均值为0、标准差为1的标准正态分布
# 之所以标准化,常用的解释是:(1)突出特征;(2)方便反向传播的计算,具体讨论可以参考:
# http://www.soolco.com/post/62169_1_1.html
T.Normalize([0.625, 0.448, 0.688],
[0.131, 0.177, 0.101]),
])

# 分块
def build_slices(self):
self.masks = []
self.files = []
self.slices = []

for i, filename in enumerate(self.csv.index.values):
filepath = (self.path /'train'/(filename+'.tiff')).as_posix()
self.files.append(filepath)

print('Transform', filename)

with rasterio.open(filepath, transform = identity) as dataset:
self.masks.append(rle_decode(self.csv.loc[filename, 'encoding'], dataset.shape))
slices = make_grid(dataset.shape, window=self.window, min_overlap=self.overlap)

# 排除掉一些没有目标对象的分块
for slc in tqdm(slices):
x1,x2,y1,y2 = slc
if self.masks[-1][x1:x2,y1:y2].sum() > self.threshold or np.random.randint(100) > 120:
self.slices.append([i,x1,x2,y1,y2])

image = dataset.read([1,2,3],
window=Window.from_slices((x1,x2),(y1,y2)))

# if image.std().mean() < 10:
# continue

# print(image.std().mean(), self.masks[-1][x1:x2,y1:y2].sum())

image = np.moveaxis(image, 0, -1)
self.x.append(image)
self.y.append(self.masks[-1][x1:x2,y1:y2])

# get data operation
def __getitem__(self, index):
image, mask = self.x[index], self.y[index]
augments = self.transform(image=image, mask=mask)
return self.as_tensor(augments['image']), augments['mask'][None]

def __len__(self):
"""
Total number of samples in the dataset
"""
return self.len

WINDOW=1024
MIN_OVERLAP=32
NEW_SIZE=256

# albumentations是一个基于OpenCV的数据增强库,拥有非常简单且强大的可以用于多种任务(分割、检测)的接口,易于定制且添加其他框架非常方便。
# 它可以对数据集进行逐像素的转换,如模糊、下采样、高斯噪声、高斯模糊、动态模糊、RGB转换、随机雾化等;
# 也可以进行空间转换(同时也会对目标标签进行转换),如裁剪、翻转、随机裁剪等
# 数据增强的方式是这样的:
# 比如在一个epoch之内,我是把所有的图片都过一遍,对于每张图片我都是进行一个transform的操作,比如transform内部有0.5的概率进行左右flip,
# 那么这张图片左右flip的概率就是0.5,可能这一个epoch不flip,下一个epoch就会flip.
# 换句话说,现有的数据增强是带有随机性的,比如是否随机镜像翻转,随机crop的时候选择哪块区域,加多强的噪声等,每次增强的结果可能都不一样,
# 这样模型相当于看到了很多份不同的图像。如果没有概率性的操作,即p=1, 做一次增强之后便不再变化,则和不做增强是等价的,即模型在整个训练过程中只能看到一份相同的不断重复的数据。
# 可以详见
# https://discuss.gluon.ai/t/topic/1666
trfm = A.Compose([
A.Resize(NEW_SIZE,NEW_SIZE),
A.HorizontalFlip(p=0.5),
A.VerticalFlip(p=0.5),

A.OneOf([
A.RandomContrast(),
A.RandomGamma(),
A.RandomBrightness(),
A.ColorJitter(brightness=0.07, contrast=0.07,
saturation=0.1, hue=0.1, always_apply=False, p=0.3),
], p=0.3),

A.OneOf([
A.ElasticTransform(alpha=120, sigma=120 * 0.05, alpha_affine=120 * 0.03),
A.GridDistortion(),
A.OpticalDistortion(distort_limit=2, shift_limit=0.5),
], p=0.0),
A.ShiftScaleRotate(),
])


# 使用模型库里的模型时,有的模型所需要的输入图像的尺寸是固定的,比如必须是64*64等,所以在图像增强后要保证图像的尺寸满足该要求。
# 可以采取以下方法来满足任意图像尺寸的输入:
#(1)传统办法:预处理,也就是边缘补0,或者切割,或者重采样、resize来得到相同大小的输入图片作为input。
#(2)模型方法:使用Kaiming He大佬的SPP-Net可以输入任意大小的图片,原文arxiv地址:https://arxiv.org/abs/1406.4729
#(3)优雅方法:加一层torch.nn.AdaptiveMaxPool2d。
# 具体讨论可以见
# https://www.zhihu.com/question/45873400

ds = HubDataset(DATA_PATH, window=WINDOW, overlap=MIN_OVERLAP, transform=trfm)

# 划分验证集和训练集
valid_idx, train_idx = [], []
for i in range(len(ds)):
# 挑出第8张图片的所有切片作为验证集
if ds.slices[i][0] == 7:
valid_idx.append(i)
else:
train_idx.append(i)

train_ds = D.Subset(ds, train_idx)
valid_ds = D.Subset(ds, valid_idx)

# define training and validation data loaders
loader = D.DataLoader(
train_ds, batch_size=BATCH_SIZE, shuffle=True, num_workers=0)

vloader = D.DataLoader(
valid_ds, batch_size=BATCH_SIZE, shuffle=False, num_workers=0)

离线制作数据集

如上所述,上述制作数据集的方式是离线进行的,该数据集准备过程参考了如下notebook:
256x256 images
具体切分方式与上面在线制作时不同,但表达的意思相同。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# 对迭代器进行tqdm封装,传入total参数来指明预计迭代次数
# https://ptorch.com/news/170.html
# 如果想选择dataframe某一行作为测试,可以参考dataframe各种索引方法
# https://www.jianshu.com/p/32bfb327bf07

for index, rles in tqdm(df_mask.iterrows(), total=len(df_mask)):
print(index)
img, mask = read_image(index)
shape = img.shape
tile_before_compress = compress * tile_size

# 计算填充量,使得可以整除
pad0 = tile_before_compress - shape[0] % tile_before_compress
pad1 = tile_before_compress - shape[1] % tile_before_compress

# 用0填充
# https://numpy.org/doc/stable/reference/generated/numpy.pad.html
# 对于img,x和y两个维度都填充,第三维度则不填充,所以第三维度设为(0, 0)
img = np.pad(img, ((pad0 // 2, pad0 - pad0 // 2), (pad1 // 2, pad1 - pad1 // 2), (0, 0)), constant_values=0)
mask = np.pad(mask, ((pad0 // 2, pad0 - pad0 // 2), (pad1 // 2, pad1 - pad1 // 2)), constant_values=0)

# 压缩图像
img = cv2.resize(img, (img.shape[1]//compress, img.shape[0]//compress), interpolation = cv2.INTER_AREA)

# 将图像分块成tile大小
img = img.reshape(img.shape[0] // tile_size, tile_size, img.shape[1] // tile_size, tile_size, 3)

# 将图像的通道调整顺序,横纵个数放在前面,然后让两者相乘,得到分块总个数
img = img.transpose(0, 2, 1, 3, 4).reshape(-1, tile_size, tile_size, 3)
mask = cv2.resize(mask, (mask.shape[1]//compress, mask.shape[0]//compress), interpolation = cv2.INTER_AREA)
mask = mask.reshape(mask.shape[0] // tile_size, tile_size, mask.shape[1] // tile_size, tile_size)
mask = mask.transpose(0, 2, 1, 3).reshape(-1, tile_size, tile_size)

# 使用zip函数将img和mask中的元素打包在一块,统一调用
# https://www.runoob.com/python3/python3-func-zip.html
for i, (im, m) in enumerate(zip(img, mask)):
x_tot.append((im / 255.0).reshape(-1, 3).mean(0))
x2_tot.append((im / 255.0).reshape(-1, 3).mean(0))
cv2.imwrite(OUT_IMG + f'{index}_{i}.png', cv2.cvtColor(im, cv2.COLOR_RGB2BGR))
cv2.imwrite(OUT_MASK + f'{index}_{i}.png', m)

img_avr = np.array(x_tot).mean(0)
img_std = np.sqrt(np.array(x2_tot).mean(0) - img_avr**2)
print("mean: ", img_avr, ", std:", img_std)

训练模型

首先通过上面的get_model()加载模型,然后再定义一系列的必要步骤,包括优化器、损失评价等,如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# 自定义Soft Dice Loss
# 网友总结了用于图像分割的常用的损失函数的PyTorch实现,见:
# https://github.com/JunMa11/SegLoss
class SoftDiceLoss(nn.Module):
def __init__(self, smooth=1., dims=(-2,-1)):
super(SoftDiceLoss, self).__init__()
self.smooth = smooth
self.dims = dims

def forward(self, x, y):
tp = (x * y).sum(self.dims)
fp = (x * (1 - y)).sum(self.dims)
fn = ((1 - x) * y).sum(self.dims)
dc = (2 * tp + self.smooth) / (2 * tp + fp + fn + self.smooth)
dc = dc.mean()
return 1 - dc

# 损失函数1
bce_fn = nn.BCEWithLogitsLoss()

# 损失函数2
dice_fn = SoftDiceLoss()

# 最终损失是两种损失函数的加权平均
def loss_fn(y_pred, y_true):
bce = bce_fn(y_pred, y_true)
dice = dice_fn(y_pred.sigmoid(), y_true)
return 0.8*bce+ 0.2*dice

# 在验证集上运行模型
@torch.no_grad()
def validation(model, loader, loss_fn):
losses = []
model.eval()

for image, target in loader:
image, target = image.to(DEVICE), target.float().to(DEVICE)
output = model(image)['out']
loss = loss_fn(output, target)
losses.append(loss.item())

return np.array(losses).mean()

### 为了让中间输出结果好看,专门定制了一个显示效果
header = r'''
Train | Valid
Epoch | Loss | Loss | Time, m
'''
# Epoch metrics time
raw_line = '{:6d}' + '\u2502{:7.3f}'*2 + '\u2502{:6.2f}'

best_loss = 10
EPOCHES = 10

# 开始迭代训练
for epoch in range(1, EPOCHES+1):
losses = []
# 计时开始
start_time = time.time()

# 这里显式地设置model是train模式,以使dropout和batchnorm等网络层中的参数不固定
# 它仅仅是一个flag,与之类似的:model.eval()是设定eval模式,即这些网络层中的参数固定,以保证测试结果的可重复性
# https://stackoverflow.com/questions/51433378/what-does-model-train-do-in-pytorch
model.train()

# 开始在训练集数据中迭代
for image, target in loader:
image, target = image.to(DEVICE), target.float().to(DEVICE)
optimizer.zero_grad()
output = model(image)['out']
loss = loss_fn(output, target)
loss.backward()
optimizer.step()
losses.append(loss.item())

# 在验证集上运行模型
vloss = validation(model, vloader, loss_fn)
print(raw_line.format(epoch, np.array(losses).mean(), vloss,
(time.time()-start_time)/60**1))
losses = []
# 及时地将最佳模型存储下来
if vloss < best_loss:
best_loss = vloss
torch.save(model.state_dict(), 'model_best.pth')

# 训练结束后删除这些存储数据地变量,并回收内存
del loader, vloader, train_ds, valid_ds, ds
gc.collect();

训练完后的模型就开始接下来在上面的测试集上进行推导,并提交结果文件(注意在kaggle上提交时,可以把前面的训练过程去掉,直接提交训练好的模型)。