【非配置】PaddleSeg DeepLabV3+实现人体肾组织图像中肾小球识别

一、项目背景

估计显示，地球上有超过70亿人，银河系中有3000亿颗恒星。相比之下，成年人体内含有37万亿个细胞。确定这些细胞之间的功能和关系是一项艰巨的任务。如果我们更好地了解细胞活动，人类健康的许多领域都会受到影响.

正如人类基因组计划绘制了整个人类DNA一样，人类生物分子图谱计划（HuBMAP）是一项重大努力。由美国国立卫生研究院（NIH）赞助，HuBMAP正在努力促进开发一个框架，该框架在历史上首次在肾小球功能组织单位水平上绘制人体图谱。HuBMAP希望成为世界上最大的合作生物项目之一，旨在成为细胞水平上人体的开放地图。
“黑客肾脏”，首先以单细胞分辨率绘制人类肾脏。

二、前期准备

安装所需的库，因为本项目是基于代码实现的，大多的PaddleSeg开源项目是使用配置文件方式实现的，但具体的API细节文档较少，所以可以提前把PaddleSeg源代码Clone下载，作参考。

1.安装PaddleSeg

!pip install paddleseg

2.安装rasterio

!pip install rasterio

3. 解压数据集

!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

三、数据集预处理

下段预处理代码来源：Kaggle

def rle_encode(im):
    pixels = im.flatten(order = 'F')
    pixels = np.concatenate([[0], pixels, [0]])
    runs = np.where(pixels[1:] != pixels[:-1])[0] + 1
    runs[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 -= 1
    ends = starts + lengths
    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)
    for lo, hi in zip(starts, ends):
        img[lo:hi] = 1
    return img.reshape(shape, order='F')

@numba.njit()
def rle_numba(pixels):
    size = len(pixels)
    points = []
    if pixels[0] == 1:
        flag = False
        points.append(1)
    else:
        flag = True
    for i in range(1, size):
        if pixels[i] != pixels[i-1]:
            if flag:
                points.append(i+1)
                flag = False
            else:
                points.append(i+1 - points[-1])
                flag = True
    if pixels[-1] == 1: points.append(size-points[-1]+1)    
    return points

def 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 = 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)

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 = 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.build_slices()
        self.len = len(self.slices)
        self.resize = T.Resize(resize)
        self.num_classes = 2
        self.img_channels = 3
        self.ignore_index = 255
        self.edge = False
    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)
            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 = slc
                    if 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 result
    
    def __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\times H\times W$的概率图，其中N是一个向量，向量的维数表示为类别的数量，也就是说向量的每个元素对应一个类别，而这个元素的数值是概率值，表示为对应该类别的概率。因为要给出每一个像素类别，所以$H\times W$大小需要与原图保持一致。分割网络可以看作是一个没有全连接层的全卷积网络，常见的网络结构可以分为Encoder和Decoder两个部分：

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

一个好的分割网络需要解决哪些问题？第一个问题就是如何通过较少Encoder层数来获得较大的感受野。首先来看下什么是感受野。输出特征图上每个点的数值，是由输入图片上大小为$k_h\times k_w$的区域的元素与卷积核每个元素相乘再相加得到的，所以输入图像上$k_h\times k_w$区域内每个元素数值的改变，都会影响输出点的像素值。我们将这个区域叫做输出特征图上对应点的感受野。感受野内每个像素数值的变动，都会影响输出点的数值变化。因此感受野也可以表示为特征图上每个点所受图像区域范围影响的大小。感受野的大小是分割网络的特征表示能力的一个重要评价指标，因为我们希望尽可能提高感受野的大小。最直接的方法就是提高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,#填写训练集的dataset
    val_dataset=valid_ds,#填写验证集的dataset
    optimizer=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.0024


2022-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的函数输出，我们可以看到以下几个评价指标。

1. mIoU: 0.9066
2. Acc: 0.9846
3. Kappa: 0.8987
4. 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 overlay


def preshow(ds,start=0):
    plt.figure(figsize=(10, 10))
    i = 0
    mask_idx = 0
    param_state_dict = paddle.load('MyOutput/best_model/model.pdparams')
    model.set_dict(param_state_dict)
    model.eval() #预测模式
    for d in ds:
        if i > 8: 
            break
        image = 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) # 创建一个和原图的一样大小的ndarray
        redimg[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 += 3
        mask_idx += 1

    plt.show()

print('从左到右，分别是原图像，标签值，预测值')
preshow(valid_ds)

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

，分别是原图像，标签值，预测值

七、总结

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

此文章为搬运
原项目链接