基于3D分割实现自动判断脂肪肝和严重程度

AI Studio

已于 2023-03-21 21:12:48 修改

阅读量321

点赞数

分类专栏：算法人工智能文章标签： 3d 人工智能机器学习

于 2023-03-14 14:57:54 首次发布

原文链接：https://aistudio.baidu.com/aistudio/projectdetail/5707698?forkThirdPart=1

版权

人工智能同时被 2 个专栏收录

180 篇文章 33 订阅

订阅专栏

算法

15 篇文章 1 订阅

订阅专栏

★★★ 本文源自AlStudio社区精品项目，【点击此处】查看更多精品内容 >>>

1.项目背景

1.1啥是脂肪肝

现在人们的日常生活方式和饮食结构发生了巨大的变化。大概就是：吃好了，动少了。

导致体内的过量的甘油三酯无法代谢，最终聚集在肝细胞内，正常肝脏就变成脂肪肝。

1.2 判断是否脂肪肝

在医学影像中，判断患者是否有脂肪肝。可以通过CT或者B超技术手段进行判断。经相关学者研究得出CT的脂肪肝检出比例高于B超，用于脂肪肝诊断时CT检出率更高，特异性更强。

在CT检查中，是通过计算肝脏与脾脏的CT值的比值来判断患者是否有脂肪肝，和脂肪肝的严重程度。参照中华医学会肝病学分会制定的标准：肝脾CT值大于1的为正常肝脏。肝脾CT比值在[0.7,1.0]间为轻度脂肪肝，在[0.5,0.7]间为中度脂肪肝，小于0.5为重度脂肪肝。

日常工作中，放射医师需要手动对肝、脾最大层面选取一定范围的ROI，然后计算两者的ROI范围内的CT总值，然后计算两者CT总值的比值，可见在CT影像上判断患者是否脂肪肝还是具有较大的工作量。

1.3如何通过深度学习实现自动判定脂肪肝和严重程度。

在日常医疗工作中，存在各种各样繁琐、重复的工作。在ct中判断脂肪肝就是一种。

现在就通过深度学习中的语义分割来辅助医生来解决手动测量脂肪肝的问题。

具体方法如下图：

2.数据准备

数据来源：医学图像十项全能挑战赛，对包含肝脏的数据，和包含脾脏的比赛数据进行整合。生产同事具有肝脏和脾脏掩膜的分割数据。

数量：88个带掩膜的数据。4个测试集。

标签：背景为0，肝脏为1，脾脏为2

格式：NIFIT格式。扩展名:nii.gz

数据用itk-snap读取展示如下：

3.MedicalSeg医疗3D分割工具

训练框架采用PaddleSeg 分割套件中的MedicalSeg医疗3D分割工具。

因为这个项目主要想用人工智能来解决医学临床问题。分割套件的使用方式不多说，具体的使用方式，可以看我这个项目

快速上手PaddleSeg的医疗3D分割套件MedicalSeg：https://aistudio.baidu.com/aistudio/projectdetail/3878920?contributionType=1&sUid=181096&shared=1&ts=1677648488476

4.分割流程

clone PaddleSeg分割套件。
解压数据。
安装分割套件对应的依赖。
对数据进行预处理。把医疗文件转换成适合模型读取的numpy文件。
开始训练
计算模型在验证集的分割得分。
导出模型方面预测。

#1.clone PaddleSeg分割套件。
!git clone https://gitee.com/PaddlePaddle/PaddleSeg.git

#2.解压数据
!unzip -o /home/aistudio/data/data194126/SpleenAndLiver.zip -d /home/aistudio/work

#3.安装分割套件对应的依赖。
%cd PaddleSeg/contrib/MedicalSeg

!pip install -r requirements.txt -i https://mirror.baidu.com/pypi/simple

#需要把自己自定义的预处理文件，复制到PaddleSeg/contrib/MedicalSeg/tools 文件中，不然无法运行
!cp /home/aistudio/prepare_SpleenAndLiver.py /home/aistudio/PaddleSeg/contrib/MedicalSeg/tools/prepare_SpleenAndLiver.py

#4. 对数据进行预处理。把医疗文件转换成适合模型读取的numpy文件。
#运行预处理文件，把SimpleITK文件转换成numpy文件，生成对应的train.txt和val.txt，和数据参数有关的json文件
!python tools/prepare_SpleenAndLiver.py

#5. 开始训练
!python3 train.py --config /home/aistudio/SpleenAndLiver.yml \
    --save_dir  "/home/aistudio/output/SpleenAndLiver_vent_128" \
    --save_interval 70 --log_iters 20 \
    --keep_checkpoint_max 4 \
    --num_workers 1 --do_eval --use_vdl

#6. 计算模型在验证集的分割得分。
#训练了3000多轮。肝脏的分割精度有0.923，脾脏的分割精度0.899
"""
2023-02-27 20:18:38 [INFO]	[EVAL] #Images: 18, Dice: 0.9391, Loss: 0.120879
2023-02-27 20:18:38 [INFO]	[EVAL] Class dice: 
[0.9946 0.9229 0.8998]
"""
!python3 val.py --config /home/aistudio/SpleenAndLiver.yml \
--model_path /home/aistudio/output/SpleenAndLiver_vent_128/best_model/model.pdparams \
--save_dir  /home/aistudio/output/SpleenAndLiver_vent_128/best_model

#7. 导出模型方面预测。
!python export.py --config /home/aistudio/SpleenAndLiver.yml \
--model_path /home/aistudio/output/SpleenAndLiver_vent_128/best_model/model.pdparams \
--save_dir /home/aistudio/export_model

5.预测，把预测结果生成NIFIT格式

MedicalSeg需要加载npy文件进行预测。生成结果也是npy文件，但是判断脂肪肝，需要在原始数据上进行。
不然会影响数据的真实性。因为需要对预测的mask结果重新采样成原始数据的大小、形状、空间等。

#预测用到的函数和类
import numpy as np 
import SimpleITK as sitk
from paddle.inference import create_predictor,Config

class Predictor:
    """
    用于预测的类
    """
    def __init__(self,model_path,param_path):
        self.pred_cfg = Config(model_path,param_path)
        self.pred_cfg.disable_glog_info()
        self.pred_cfg.enable_memory_optim()
        self.pred_cfg.switch_ir_optim(True)
        self.pred_cfg.enable_use_gpu(100, 0)
        # self.pred_cfg.disable_gpu()
        self.predictor = create_predictor(self.pred_cfg)

    def predict(self, data):
        input_names = self.predictor.get_input_names()
        input_handle = self.predictor.get_input_handle(input_names[0])
        output_names = self.predictor.get_output_names()
        output_handle = self.predictor.get_output_handle(output_names[0])
        input_handle.reshape(data.shape)
        input_handle.copy_from_cpu(data)
        self.predictor.run()
        result = output_handle.copy_to_cpu()
        return result

def resampleImage(sitkimg,new_shape,new_spacing):
    #对SimpleITK 的数据进行重新采样。重新设置spacing和shape
    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(sitkimg)  
    resampler.SetOutputSpacing(new_spacing)
    resampler.SetSize(new_shape)
    resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
    resampler.SetInterpolator(sitk.sitkLinear)
    return resampler.Execute(sitkimg)  

def crop_wwwc(sitkimg,max_v,min_v):
    #对SimpleITK的数据进行窗宽窗位的裁剪，应与训练前对数据预处理时一致
    intensityWindow = sitk.IntensityWindowingImageFilter()
    intensityWindow.SetWindowMaximum(max_v)
    intensityWindow.SetWindowMinimum(min_v)
    return intensityWindow.Execute(sitkimg)

def GetLargestConnectedCompont(binarysitk_image):
    # 最大连通域提取,binarysitk_image 是掩膜
    cc = sitk.ConnectedComponent(binarysitk_image)
    stats = sitk.LabelIntensityStatisticsImageFilter()
    stats.SetGlobalDefaultNumberOfThreads(8)
    stats.Execute(cc, binarysitk_image)#根据掩膜计算统计量
    # stats.
    maxlabel = 0
    maxsize = 0
    for l in stats.GetLabels():#掩膜中存在的标签类别
        size = stats.GetPhysicalSize(l)
        if maxsize < size:#只保留最大的标签类别
            maxlabel = l
            maxsize = size
    labelmaskimage = sitk.GetArrayFromImage(cc)
    outmask = labelmaskimage.copy()
    if len(stats.GetLabels()):
        outmask[labelmaskimage == maxlabel] = 255
        outmask[labelmaskimage != maxlabel] = 0
    return outmask

#GPU下进行预测

origin_path = '/home/aistudio/work/SpleenAndLiver/test/liver_91_0000.nii.gz'
origin = sitk.ReadImage(origin_path)

new_shape = (128, 128, 128) #xyz #这个形状与训练的对数据预处理的形状要一致
image_shape = origin.GetSize()
spacing = origin.GetSpacing()
new_spacing = tuple((image_shape / np.array(new_shape)) *spacing) 

itk_img_res = resampleImage(origin,new_shape,new_spacing)  # 得到重新采样后的图像
itk_img_res = crop_wwwc(itk_img_res,max_v=300,min_v=-100)#和预处理文件一致
npy_img = sitk.GetArrayFromImage(itk_img_res).astype("float32")
input_data = np.expand_dims(npy_img,axis=0)
if input_data.max() > 0: #归一化
    input_data = input_data / input_data.max()
input_data = np.expand_dims(input_data,axis=0)
print(f"输入网络前数据的形状:{input_data.shape}")#shape(1, 1, 128, 128, 256)

#创建预测器，加载模型进行预测
predictor = Predictor('/home/aistudio/export_model/model.pdmodel',
                        '/home/aistudio/export_model/model.pdiparams')
output_data = predictor.predict(input_data)
print(f"预测结果的形状：{output_data.shape}")#shape (1, 128, 128, 256)

#加载3d模型预测的mask，由numpy 转换成SimpleITK格式
data = np.squeeze(output_data)
mask_itk_new = sitk.GetImageFromArray(data)
mask_itk_new.SetSpacing(new_spacing)
mask_itk_new.SetOrigin(origin.GetOrigin())
mask_itk_new.SetDirection(origin.GetDirection())
mask_itk_new = sitk.Cast(mask_itk_new,sitk.sitkUInt8)

x,y,z = mask_itk_new.GetSize()
mask_array = np.zeros((z,y,x),np.uint8)
max_value = np.max(sitk.GetArrayViewFromImage(mask_itk_new))
#对转换成SimpleITK的预测mask进行处理，只保留最大连通域，去除小目标
for index in range(1,max_value+1):
    sitk_seg = sitk.BinaryThreshold(mask_itk_new, lowerThreshold=index, upperThreshold=index, insideValue=255, outsideValue=0)
    # step2.形态学开运算
    BMO = sitk.BinaryMorphologicalOpeningImageFilter()
    BMO.SetKernelType(sitk.sitkNearestNeighbor)
    BMO.SetKernelRadius(2)
    BMO.SetForegroundValue(1)
    sitk_open = BMO.Execute(sitk_seg!=0)
    #提取每个椎体的最大连通域提取，为了去掉小目标
    sitk_open_array = GetLargestConnectedCompont(sitk_open)
    mask_array[sitk_open_array==255] = int(index)

#对处理好的预测mask，重采样原始的size 和spacing
sitkMask = sitk.GetImageFromArray(mask_array)
sitkMask.CopyInformation(mask_itk_new)
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(sitkMask)  # 需要重新采样的目标图像
resampler.SetSize(origin.GetSize())
resampler.SetOutputSpacing(origin.GetSpacing())
resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
resampler.SetInterpolator(sitk.sitkNearestNeighbor)

Mask = resampler.Execute(sitkMask)  # 得到重新采样后的图像
Mask.CopyInformation(origin)
sitk.WriteImage(Mask,'/home/aistudio/pred_data.nii.gz')
print("预测成功！")

使用itk-snap加载标签文件和预测数据文件，看到模型有较好的预测能力，但是细节上达不到很好的水平。
可以增加样本、更换模型，增加训练轮次来改善模型精度。

6. 自动判断是否脂肪肝和脂肪肝的严重程度

上面已经对原始数据进行预测，并生成mask，然后再转换成原始数据同样的参数。

现在有处理好的mask文件了。

可以开始对肝脏和脾脏随机割立方体，来计算平均比值。

具体方法如下：

根据mask结果，缩小肝脏和脾脏的范围。
先对肝脏随机获取5个立方体。
再对脾脏随机获取5个立方体。
5个肝脏立方体和5个脾脏立方体。两两配对成30多对子。
计算每个对子的CT比值。再计算所有CT比值的平均CT比值。
根据医学标准来划分是否有脂肪肝，和脂肪肝的严重程度。

import SimpleITK as sitk
import random
import numpy as np
import itertools
def maskcroppingbox(mask):
    #寻找mask范围的3D最大边界范围
    mask_2 = np.argwhere(mask)
    (zstart, ystart, xstart), (zstop, ystop, xstop) = mask_2.min(axis=0), mask_2.max(axis=0) + 1
    #让变量好看点，保证max是最大值
    if zstart- zstop < 0 :  zmax,zmin = zstop,zstart
    if ystart- ystop < 0 :  ymax,ymin = ystop,ystart
    if xstart- xstop < 0 :  xmax,xmin = xstop,xstart
    return zmax,zmin,ymax,ymin,xmax,xmin

#读取原始数据 和预测的mask数据
imgSitk = sitk.ReadImage(origin_path)
maskSitk = sitk.ReadImage('/home/aistudio/pred_data.nii.gz')

imgNp = sitk.GetArrayFromImage(imgSitk)
maskNp = sitk.GetArrayFromImage(maskSitk)

#先定义需要扣出 变成为10mm的立方体
#因为CT数据每个体素之间都有物理空间的。以mm为单位。
length = 10 

#label是一个字典。键是类别标签。值是一个列表，用来放对应类别的立方体的ct值总数
#例如要肝脏和脾脏各随机割5个立方体。label={1:[100,200,300,200,299],2:[300,200,300,200,110]}
label  = np.unique(maskNp)
label = sorted(label, reverse=False)
label.remove(0)
label = {1:list(),2:list()}

#循环肝脏 和脾脏，割完肝脏，再割脾脏
for index,value in label.items():
    print(f"开始处理类别：{index}")
    while len(value) <6: #设置6，代表我打算割五个立方体
        temp = maskNp.copy()
        temp[temp!=index] = 0
        temp[temp>0] = 1

        spacingx,spacingy,spacingz = maskSitk.GetSpacing()
        # spacing物理单位是mm，因此立方体的长度除以spacing就可得需要多少个像素点。
        lengthx = round(length/spacingx)
        lengthy = round(length/spacingy)
        lengthz = round(length/spacingz)

        #找到器官的掩膜的最大边界范围。缩小范围
        zmax,zmin,ymax,ymin,xmax,xmin = maskcroppingbox(temp)

        #在3Dmask中随机获取个坐标。
        posZ = random.randint(0,zmax-zmin)+zmin
        posY = random.randint(0,ymax-ymin)+ymin
        posX = random.randint(0,xmax-xmin)+xmin
        #这个坐标是像素点坐标。坐标为中心，生成立方体的坐标，从mask中割出来
        sliceZ = slice(posZ,posZ+lengthz)
        sliceY = slice(posY,posY+lengthy)
        sliceX = slice(posX,posX+lengthx)
        block = temp[sliceZ,sliceY,sliceX]

        #如果这个立方体都在3Dmask里面，那np.sum（）等于lengthx*lengthy*lengthz。因为我设置了mask的值都是1
        if np.sum(block) == lengthx*lengthy*lengthz:
            value.append(np.sum(imgNp[sliceZ,sliceY,sliceX]))
            print(f"{len(value)}个立方体的CT总值:{value}")

#现在有5个肝脏的立方体和5个脾脏的立方体，找出两个的所有组合，然后计算它们比值
tupleNums = list(itertools.product(label[1], label[2]))
cts = [tupleNum[0]/tupleNum[1] for tupleNum in tupleNums]
#再对这个比较做平均值
mean_ct = np.mean(cts)
print("\n*********************************")
print(f"平均CT比值为{mean_ct}")
if mean_ct <1:
    if mean_ct <1.0 and  mean_ct >=0.7:
        print("结果为:轻度脂肪肝")
    elif mean_ct <0.7 and mean_ct >= 0.5:
        print("结果为:中度脂肪肝")
    else:
        print("结果为:重度脂肪肝")
else:
    print("结果为:非脂肪肝")

开始处理类别：1
1个立方体的CT总值:[205964.0]
2个立方体的CT总值:[205964.0, 224356.0]
3个立方体的CT总值:[205964.0, 224356.0, 212064.0]
4个立方体的CT总值:[205964.0, 224356.0, 212064.0, 209330.0]
5个立方体的CT总值:[205964.0, 224356.0, 212064.0, 209330.0, 202258.0]
6个立方体的CT总值:[205964.0, 224356.0, 212064.0, 209330.0, 202258.0, 209558.0]
开始处理类别：2
1个立方体的CT总值:[203052.0]
2个立方体的CT总值:[203052.0, 209845.0]
3个立方体的CT总值:[203052.0, 209845.0, 221099.0]
4个立方体的CT总值:[203052.0, 209845.0, 221099.0, 221029.0]
5个立方体的CT总值:[203052.0, 209845.0, 221099.0, 221029.0, 223527.0]
6个立方体的CT总值:[203052.0, 209845.0, 221099.0, 221029.0, 223527.0, 206339.0]

*********************************
平均CT比值为0.9847654700279236
结果为:轻度脂肪肝
6个立方体的CT总值:[203052.0, 209845.0, 221099.0, 221029.0, 223527.0, 206339.0]

*********************************
平均CT比值为0.9847654700279236
结果为:轻度脂肪肝