基于飞桨 DeepLabV3+实现人体肾组织图像中肾小球识别

一、项目背景

估计显示，地球上有超过70亿人，银河系中有3000亿颗恒星。相比之下，成年人体内含有37万亿个细胞。确定这些细胞之间的功能和关系是一项艰巨的任务。如果我们更好地了解细胞活动，人类健康的许多领域都会受到影响.
正如人类基因组计划绘制了整个人类DNA一样，人类生物分子图谱计划（HuBMAP）是一项重大努力。由美国国立卫生研究院（NIH）赞助，HuBMAP正在努力促进开发一个框架，该框架在历史上首次在肾小球功能组织单位水平上绘制人体图谱。HuBMAP希望成为世界上最大的合作生物项目之一，旨在成为细胞水平上人体的开放地图。
“黑客肾脏”，首先以单细胞分辨率绘制人类肾脏。

二、前期准备

安装所需的库，本项目是在开源代码上实现的，所以可以提前把源代码Clone下载，这是一个非常强大的开源库，内置45+模型算法及140+预训练模型，支持配置化驱动和API调用开发方式，打通数据标注、模型开发、训练、压缩、部署的全流程，提供语义分割、交互式分割、Matting、全景分割四大分割能力，助力算法在医疗、工业、遥感、娱乐等场景落地应用，十分值得学习参考。

1.安装PaddleSeg

!pip install paddleseg

2.安装rasterio

!pip install rasterio

3. 数据集

数据集来源：HuBMAP - Hacking the Kidney
这个数据集来源于Kaggle新赛 | HuBMAP: 识别人体肾脏组织图像中的肾小球~总奖金 6 万美金
数据集高速下载链接：https://aistudio.baidu.com/aistudio/datasetdetail/178466
这是关于入侵肾脏的数据集，它的标签是二个类别，我们先把解压：

!unzip -oq /home/aistudio/data/data178466/入侵肾脏.zip -d /home/aistudio/work/dataset/

4.导入相关库

from PIL import Image as PilImage
import paddle
import numpy as np
from paddle.io import Dataset
from paddle.vision.transforms import transforms as T
import scipy.misc as misc
import pandas as pd
import pathlib, sys, os, random,io,time
import numba, cv2, gc
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
import rasterio
from rasterio.windows import Window
import paddleseg.transforms as segT
from paddleseg.datasets import Dataset
from paddleseg.core import train
from paddleseg.core import evaluate
from paddleseg.core import predict
from paddleseg.models import DeepLabV3P
from paddleseg.models import GCNet
from paddleseg.models import FCN
from paddleseg.models.backbones import HRNet_W48
from paddleseg.models.backbones import resnet_vd
from paddleseg.models.losses import CrossEntropyLoss

三、数据集预处理

def rle_encode(im):pixels = im.flatten(order = 'F')pixels = np.concatenate([[0], pixels, [0]])runs = np.where(pixels[1:] != pixels[:-1])[0] + 1runs[1::2] -= runs[::2]return ' '.join(str(x) for x in runs)def rle_decode(mask_rle, shape=(256, 256)):s = mask_rle.split()starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]starts -= 1ends = starts + lengthsimg = np.zeros(shape[0]*shape[1], dtype=np.uint8)for lo, hi in zip(starts, ends):img[lo:hi] = 1return img.reshape(shape, order='F')@numba.njit()
def rle_numba(pixels):size = len(pixels)points = []if pixels[0] == 1:flag = Falsepoints.append(1)else:flag = Truefor i in range(1, size):if pixels[i] != pixels[i-1]:if flag:points.append(i+1)flag = Falseelse:points.append(i+1 - points[-1])flag = Trueif pixels[-1] == 1: points.append(size-points[-1]+1)    return pointsdef rle_numba_encode(image):pixels = image.flatten(order = 'F')points = rle_numba(pixels)return ' '.join(str(x) for x in points)def make_grid(shape, window=256, min_overlap=32):x, y = shapenx = x // (window - min_overlap) + 1x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)x1[-1] = x - windowx2 = (x1 + window).clip(0, x)ny = y // (window - min_overlap) + 1y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)y1[-1] = y - windowy2 = (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)

1.定义Dataset

在这里，使用了自定义Dataset，因为源数据集是非常大的图片，所以每次getitem一次，都会重新读取源文件再resize成想要的大小，这是非常耗时的操作，模型训练的时候，这块也会严重影响训练时间，所以可以在这里优化。

在转换后，把转换好的数据缓存，下次需要使用直接读取缓存，将会节约很多的时间。

identity = rasterio.Affine(1, 0, 0, 0, 1, 0)class HubDataset(Dataset):def __init__(self, root_dir, transform,window=256, overlap=32, threshold = 100,resize = (256,256)):self.path = pathlib.Path(root_dir)self.overlap = overlapself.window = windowself.transform = transformself.csv = pd.read_csv((self.path / 'train.csv').as_posix(),index_col=[0])self.threshold = thresholdself.build_slices()self.len = len(self.slices)self.resize = T.Resize(resize)self.num_classes = 2self.img_channels = 3self.ignore_index = 255self.edge = Falsedef 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)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 slices:x1,x2,y1,y2 = slcif self.masks[-1][x1:x2,y1:y2].sum() > self.threshold:self.slices.append([i,x1,x2,y1,y2])def __getitem__(self, index):exists = os.path.exists(f'work/dataseta/{index}.npy')if not exists:idx = self.slices[index][0]filename = self.files[idx]x1,x2,y1,y2 = self.slices[index][1:] with rasterio.open(filename, transform = identity) as dataset:image = dataset.read([1,2,3],window=Window.from_slices((x1,x2),(y1,y2)))image = np.moveaxis(image, 0, -1)mask = self.masks[idx][x1:x2,y1:y2]augments = {"image":image, "mask":mask}result = {}result['img'] =  self.resize(np.array(augments['image']))result['label'] = self.resize(augments['mask'][None].transpose(1,2,0)).squeeze(axis=2) if not os.path.exists('work/dataseta/'):os.mkdir('work/dataseta/')np.save(f'work/dataseta/{index}.npy',result)else:result = np.load(f'work/dataseta/{index}.npy',allow_pickle=True).item()result = self.transform(result)return resultdef __len__(self):return self.len

2.创建Dataset

WINDOW=1024
MIN_OVERLAP=32
DATA_PATH = '/home/aistudio/work/dataset'
trfm = segT.Compose([segT.RandomPaddingCrop(crop_size=(256,256)),segT.RandomDistort(brightness_range=0.4,contrast_range=0.4,saturation_range=0.4),segT.Normalize()
])
ds = HubDataset(DATA_PATH, window=WINDOW, overlap=MIN_OVERLAP, transform=trfm)

3.数据集抽样展示

idx = 4
image = ds[idx]['img']
mask = ds[idx]['label']
plt.figure(figsize=(8,16))
plt.subplot(211)
plt.imshow(mask, cmap='gray')
plt.subplot(212)
plt.imshow(image[0])
_ = rle_numba_encode(mask[0])

/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/rasterio/__init__.py:220: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix be returned.s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)

4.划分数据集及训练集

生成一组序列数组，随机选择90%的序列，作训练的序列组；再用原序列差集运算训练序列，作验证集序列。

valid_idx, train_idx = [], []
trancount = (int)(len(ds)*0.9)
allindex = np.arange(len(ds))
trainindex = np.random.choice(allindex,trancount,replace=False)
valindex = np.setdiff1d(allindex,trainindex)

构造一个数据集级的数据子数据集。

给定原数据集合的指标数组，可以以此数组构造原数据集合的子数据集合。

train_ds = paddle.io.Subset(ds,trainindex)
valid_ds = paddle.io.Subset(ds, valindex)

四、DeepLabV3+模型

1.原理介绍

对于当前十分流行的实时分割场景，例如短视频领域、直播领域，都需要实时分割算法，这就需要我们的算法足够轻量。那么一个好的分割模型网络是什么样子的呢？首先我们来看下常见的分割网络是什么样的结构。
一个分割网络的输出是一个 N × H × W N\times H\times W N×H×W的概率图，其中N是一个向量，向量的维数表示为类别的数量，也就是说向量的每个元素对应一个类别，而这个元素的数值是概率值，表示为对应该类别的概率。因为要给出每一个像素类别，所以 H × W H\times W H×W大小需要与原图保持一致。分割网络可以看作是一个没有全连接层的全卷积网络，常见的网络结构可以分为Encoder和Decoder两个部分：

Encoder部分特征图的长宽会逐步缩小，也就是需要下采样。
Decoder部分特征图的长宽会逐步增大，也就是需要上采样。

一个好的分割网络需要解决哪些问题？第一个问题就是如何通过较少Encoder层数来获得较大的感受野。首先来看下什么是感受野。输出特征图上每个点的数值，是由输入图片上大小为 k h × k w k_h\times k_w kh×kw的区域的元素与卷积核每个元素相乘再相加得到的，所以输入图像上 k h × k w k_h\times k_w kh×kw区域内每个元素数值的改变，都会影响输出点的像素值。我们将这个区域叫做输出特征图上对应点的感受野。感受野内每个像素数值的变动，都会影响输出点的数值变化。因此感受野也可以表示为特征图上每个点所受图像区域范围影响的大小。感受野的大小是分割网络的特征表示能力的一个重要评价指标，因为我们希望尽可能提高感受野的大小。最直接的方法就是提高Encoder网络的层数，但是这种方法会在扩大感受野的同时，带来一些弊端，例如增加了网络参数规模，导致训练耗时增加。那么有没有更好的方法呢？

通过扩大卷积核可触达的区域来增大感受野，扩大卷积所包含的信息范围。如下图所示，该图为空洞卷积作用的示意图示意图。可以看到未使用空洞卷积的网络，一个输出的深蓝色的点可以覆盖5个输入像素点。在使用空洞卷积之后，一个输出点可以覆盖9个输入像素的点。其具体应用示例可以参见下面DeepLabv3+的网络结构。

分割网络需要面临的第二个问题是下采样过程中特征的长和宽在不断减少，我们虽然得到了高阶特征，但是维度较小，如果直接作为上采样输入，则有可能丢失一部分信息，怎么解决这个问题呢？这时可以使用跳跃连接功能。通过跳跃连接，使上采样过程不仅使用高阶特征，也可以使用同维度的低阶特征，具体使用示例可以参见基本原理中的U-Net模型介绍。

除了上述两个问题之外，分割网络还需要解决多尺度问题，多尺度问题是指同一个图片中目标有大有小，如何保证不同尺度的目标均能被很好的处理呢？这里可以使用ASPP网络结构，通俗的讲，其实现思想是让不同尺度的目标交给不同感受野的卷积层来处理，大感受野处理大目标，小感受野处理小目标。其具体应用示例可以参见下面DeepLabv3+的网络结构。

2.DeepLabv3+

DeepLabv3+是DeepLab系列的最后一篇文章，其前作有DeepLabv1、DeepLabv2和DeepLabv3。在最新作中，作者结合编码器-解码器(encoder-decoder)结构和空间金字塔池化模块(Spatial Pyramid Pooling, SPP)的优点提出新的语义分割网络DeepLabv3+，在 PASCAL VOC 2012和Cityscapes数据集上取得新的state-of-art performance.
其整体结构如下所示，Encoder的主体是带有空洞卷积(Atrous Convolution)的骨干网络，骨干网络可采用ResNet等常用的分类网络，作者使用了改进的Xception模型作为骨干网络。紧跟其后的空洞空间金字塔池化模块(Atrous Spatial Pyramid Pooling, ASPP)则引入了多尺度信息。相比前作DeepLabv3，DeepLabv3+加入decoder模块，将浅层特征和深层特征进一步融合，优化分割效果，尤其是目标边缘的效果。此外，作者将深度可分离卷积(Depthwise Separable Convolution)应用到ASPP和Decoder模块，提高了语义分割的健壮性和运行速率。

图1 DeepLabv3+结构图

具体原理细节请参考Encoder-Decoder with Atrous Separable Convolution for Semantic Image Segmentation。

net = resnet_vd.ResNet50_vd(multi_grid=[1,2,4],pretrained='https://bj.bcebos.com/paddleseg/dygraph/resnet50_vd_ssld_v2.tar.gz')
model = DeepLabV3P(2,net,aspp_ratios=[1,12,24,36])
# model = FCN(2,HRNet_W48(pretrained='https://bj.bcebos.com/paddleseg/dygraph/hrnet_w48_ssld.tar.gz')) 也可以用这个模型

W1216 11:39:10.258989  5883 gpu_resources.cc:61] Please NOTE: device: 0, GPU Compute Capability: 7.0, Driver API Version: 11.2, Runtime API Version: 11.2
W1216 11:39:10.262938  5883 gpu_resources.cc:91] device: 0, cuDNN Version: 8.2.2022-12-16 11:39:11 [INFO]   Loading pretrained model from https://bj.bcebos.com/paddleseg/dygraph/resnet50_vd_ssld_v2.tar.gz
2022-12-16 11:39:12 [INFO]  There are 275/275 variables loaded into ResNet_vd.

五、模型训练

# 设置学习率
base_lr = 0.01
lr = paddle.optimizer.lr.PolynomialDecay(base_lr,decay_steps=20000, power=0.9, end_lr=0)
optimizer= paddle.optimizer.Momentum(learning_rate=lr,parameters=model.parameters(),momentum=0.9,weight_decay=4.0e-5)
losses = {}
losses['types'] = [CrossEntropyLoss()]
losses['coef'] = [1]

因为PaddleSeg中Dataset基类已经定义了以下变量，该变量在后续的评测中会有用到，而我们目前使用的Dataset paddle.io.Subset类并没有定义，所以在这里我们定义下以下变量。

valid_ds.num_classes = 2
valid_ds.img_channels = 3
valid_ds.ignore_index = 255
valid_ds.edge = False
train_ds.num_classes = 2

train(model=model,train_dataset=train_ds,#填写训练集的datasetval_dataset=valid_ds,#填写验证集的datasetoptimizer=optimizer,#优化器save_dir='/home/aistudio/MyOutput',#保存路径batch_size=20,iters=20000,#训练次数save_interval=500,#保存的间隔次数log_iters=50,#日志打印间隔num_workers=4,losses=losses,#传入 loss函数use_vdl=True)#是否使用visualDL

六、模型验证及可视化预测结果

evaluate(model, valid_ds, num_workers=2)

2022-12-16 15:28:06 [INFO]   Start evaluating (total_samples: 177, total_iters: 177)...177/177 [==============================] - 4s 23ms/step - batch_cost: 0.0228 - reader cost: 0.00242022-12-16 15:28:11 [INFO]    [EVAL] #Images: 177 mIoU: 0.9066 Acc: 0.9846 Kappa: 0.8987 Dice: 0.9493
2022-12-16 15:28:11 [INFO]  [EVAL] Class IoU:
[0.9833 0.8299]
2022-12-16 15:28:11 [INFO]  [EVAL] Class Precision:
[0.9903 0.9202]
2022-12-16 15:28:11 [INFO]  [EVAL] Class Recall:
[0.9929 0.8943](0.9066375237795998,0.9845937093098959,array([0.98334124, 0.8299338 ]),array([0.99032686, 0.92015978]),0.8986668062498686)

从evaluate的函数输出，我们可以看到以下几个评价指标。

mIoU: 0.9066
Acc: 0.9846
Kappa: 0.8987
Dice: 0.9493

看着还不错，接下来我们再打印一下，看看模型的预测结果与标签对比。

def get_overlay(image, colored_mask):overlay = cv2.addWeighted(image.transpose(1,2,0), 1, colored_mask.transpose(1,2,0), 0.3, 0)return overlaydef preshow(ds,start=0):plt.figure(figsize=(10, 10))i = 0mask_idx = 0param_state_dict = paddle.load('MyOutput/best_model/model.pdparams')model.set_dict(param_state_dict)model.eval() #预测模式for d in ds:if i > 8: breakimage = d['img']label = d['label']cimg = np.expand_dims(image,axis=0)data = model(paddle.to_tensor(cimg))[0][0]mask = np.argmax(data, axis=0)redimg = np.zeros(image.shape,dtype=np.uint8) # 创建一个和原图的一样大小的ndarrayredimg[0,:,:][mask==1]=255 plt.subplot(3, 3, i + 1)plt.imshow(image[0])plt.axis("off")plt.subplot(3, 3, i + 2)plt.imshow(label)plt.axis("off")plt.subplot(3, 3, i + 3)plt.imshow(redimg[0])plt.axis("off")i += 3mask_idx += 1plt.show()print('从左到右，分别是原图像，标签值，预测值')
preshow(valid_ds)

从左到右，分别是原图像，标签值，预测值

七、总结

本项目的完成，花了不少时间，一开始我使用自己的方式完成了这个项目，但评测mIou只有60%，我又尝试PaddleSeg配置文件方式，发现mIou有90%，我就纳闷为什么会相差这么多，于是我Debuge PaddleSeg一步一步还原deeplabv3p_resnet50_os8_voc12aug_512x512_40k.yml这个配置文件的Dataset、Transforms、学习率和优化器等，完成了此项目，收获很大，更加深入的了解PaddleSeg的底层。