一. 环境
- 安装基本环境(自行安装)
- 创建python的虚拟环境:
conda create -n nnUnet python=3.9
- 查看服务器版本,开GPU
nvidia-smi
- 安装对应版本的Pytorch:
-
官网:Pytorch版本
-
复制、运行
conda install pytorch torchvision torchaudio pytorch-cuda=11.8 -c pytorch -c nvidia
- 查看服务器环境
-
第一种方法:运行程序
-
第二种方法:直接查看
import torchvision
import torch
print(torch.version.cuda)
print(torch.__version__)
print(torchvision.__version__)
print(torch.cuda.is_available())
# 11.3
# 1.12.0
# 0.13.0
# True
二. 安装nnunet:两种方法
- 下载nnunet。两种方法任选一
git clone git://github.com/MIC-DKFZ/nnUNet.git
- 下载nnunet.zip。百度云提取码:f8xf
- 进入nnUNet文件夹
cd nnUNet
- 安装所需依赖包,向终端添加新的命令,这些命令用于后续整个nnU-Net pipeline的执行,这些命令都有一个前缀:nnUNetv2_ (两种方法任选一)
- 第一种方法:
pip install -e .
最好用镜像安装,速度快,如下代码:
pip install -e . -i https://pypi.tuna.tsinghua.edu.cn/simple
- 第二种方法:pip install -e .相当于python setup.py,也就是运行文件夹中的setup.py文件
python setup.py install
- 安装隐藏层(可选,建议安装),hiddenlayer 使 nnU-net 能够生成网络拓扑图。
- 第一种方法:终端运行
pip install --upgrade git+https://github.com/FabianIsensee/hiddenlayer.git
- 第二种方法:
- 下载hiddenlayer.zip,官网链接https://github.com/FabianIsensee/hiddenlayer.git,或是直接下载的百度云提取码:ysoq
- 解压zip,查看setup.py
- 安装hiddenlayer
python setup.py install
三. 数据集结构化处理
3.1 新建文件夹
在nnUNet-wh文件夹下,新建DATASET。
- 依次新建文件夹:dataset_conversion、nnUNet_preprocessed、nnUNet_raw、nnUNet_trained_models。
- 在nnUNet_raw文件中,新建文件夹Dataset040_KiTS。
3.2 设置 nnUNet 读取文件的路径
- 命令行,确保在nnunet激活环境下,输入vim ~/.bashrc,然后点击键盘insert开始插入,在bashrc文末,添加如下三行代码。然后按住Esc键,输入:再输入wq,就可以保存退出。
vim ~/.bashrc
'''
说明,这里是路径是你自己的路径,就是上一步创建的三个文件夹的路径(这部分说明不需要写进去,只需要以下三行代码)
'''
export nnUNet_raw="/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_raw"
export nnUNet_preprocessed="/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_preprocessed"
export nnUNet_results="/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_trained_models"
然后命令行输入source ~/.bashrc,确保激活路径。
重点: 然后分别键入三个echo $nnUNet_results,验证是否可以识别。不能识别,后续无法进行数据预处理。
echo $nnUNet_results
echo $nnUNet_raw
echo $nnUNet_preprocessed
3.3 数据集重命名
- 在dataset_conversion文件中,新建一个Dataset040_KiTS.py文件
第一种方法代码:
import os
import json
import shutil
def save_json(obj, file, indent=4, sort_keys=True):
with open(file, 'w') as f:
json.dump(obj, f, sort_keys=sort_keys, indent=indent)
def maybe_mkdir_p(directory):
directory = os.path.abspath(directory)
splits = directory.split("/")[1:]
for i in range(0, len(splits)):
if not os.path.isdir(os.path.join("/", *splits[:i + 1])):
try:
os.mkdir(os.path.join("/", *splits[:i + 1]))
except FileExistsError:
# this can sometimes happen when two jobs try to create the same directory at the same time,
# especially on network drives.
print("WARNING: Folder %s already existed and does not need to be created" % directory)
def subdirs(folder, join=True, prefix=None, suffix=None, sort=True):
if join:
l = os.path.join
else:
l = lambda x, y: y
res = [l(folder, i) for i in os.listdir(folder) if os.path.isdir(os.path.join(folder, i))
and (prefix is None or i.startswith(prefix))
and (suffix is None or i.endswith(suffix))]
if sort:
res.sort()
return res
base = "/root/autodl-tmp/kits" # 原始数据集路径
out = "/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_raw/Dataset040_KiTS" # 结构化数据集目录
cases = subdirs(base, join=False)
maybe_mkdir_p(out)
maybe_mkdir_p(os.path.join(out, "imagesTr"))
maybe_mkdir_p(os.path.join(out, "imagesTs"))
maybe_mkdir_p(os.path.join(out, "labelsTr"))
for c in cases:
case_id = int(c.split("_")[-1])
if case_id < 210:
shutil.copy(os.path.join(base, c, "imaging.nii.gz"), os.path.join(out, "imagesTr", c + "_0000.nii.gz"))
shutil.copy(os.path.join(base, c, "segmentation.nii.gz"), os.path.join(out, "labelsTr", c + ".nii.gz"))
else:
shutil.copy(os.path.join(base, c, "imaging.nii.gz"), os.path.join(out, "imagesTs", c + "_0000.nii.gz"))
json_dict = {}
"""
name: 数据集名字
dexcription: 对数据集的描述
modality: 模态,0表示CT数据,1表示MR数据。nnU-Net会根据不同模态进行不同的预处理(nnunet-v2版本改为channel_names)
labels: label中,不同的数值代表的类别(v1版本和v2版本的键值对刚好是反过来的)
file_ending: nnunet v2新加的
numTraining: 训练集数量
numTest: 测试集数量
training: 训练集的image 和 label 地址对
test: 只包含测试集的image. 这里跟Training不一样
"""
json_dict['name'] = "KiTS"
json_dict['description'] = "kidney and kidney tumor segmentation"
json_dict['tensorImageSize'] = "4D"
json_dict['reference'] = "KiTS data for nnunet"
json_dict['licence'] = ""
json_dict['release'] = "0.0"
json_dict['channel_names'] = {
"0": "CT",
}
json_dict['labels'] = {
"background": "0",
"Kidney": "1",
"Tumor": "2"
}
json_dict['numTraining'] = len(cases) # 应该是210例,直接写210
json_dict['file_ending'] = ".nii.gz"
json_dict['numTest'] = 0
json_dict['training'] = [{'image': "./imagesTr/%s.nii.gz" % i, "label": "./labelsTr/%s.nii.gz" % i} for i in cases]
#json_dict['test'] = []
save_json(json_dict, os.path.join(out, "dataset.json"))
这个方法生成的文件如图:
第二种方法代码:
我在autodl-fs/nnUNet-wh/DATASET/dataset_conversion/Dataset210_KiTS2019.py
路径下新建了一个py文件,复制nnUNetV2版本中的autodl-fs/nnUNet-wh/nnunetv2/dataset_conversion/Dataset220_KiTS2023.py
,需要修改标签和数量。具体修改内容如下:
将标签修改为[0 1 2],因为kits19只有背景0,肾脏1和癌症2,用于训练的数量为210个,
而kits23多了一个标签3,用于训练的数据为220个。
from batchgenerators.utilities.file_and_folder_operations import *
import shutil
from nnunetv2.dataset_conversion.generate_dataset_json import generate_dataset_json
from nnunetv2.paths import nnUNet_raw
def convert_kits2023(kits_base_dir: str, nnunet_dataset_id: int = 209):
task_name = "KiTS2019"
foldername = "Dataset%03.0d_%s" % (nnunet_dataset_id, task_name) # 生成的文件名字:Dataset209_KiTS2019
# setting up nnU-Net folders
out_base = join(nnUNet_raw, foldername)
imagestr = join(out_base, "imagesTr")
labelstr = join(out_base, "labelsTr")
# imagests = join(out_base, "imagesTs") # 生成测试集
maybe_mkdir_p(imagestr)
maybe_mkdir_p(labelstr)
maybe_mkdir_p(imagests)
cases = subdirs(kits_base_dir, prefix='case_', join=False)
for tr in cases:
case_id = int(tr.split("_")[-1])
if case_id < 210:
shutil.copy(join(kits_base_dir, tr, 'imaging.nii.gz'), join(imagestr, f'{tr}_0000.nii.gz'))
shutil.copy(join(kits_base_dir, tr, 'segmentation.nii.gz'), join(labelstr, f'{tr}.nii.gz'))
else:
pass
# shutil.copy(join(kits_base_dir, tr, 'imaging.nii.gz'), join(imagests, f'{tr}_0000.nii.gz')) #
generate_dataset_json(out_base, {0: "CT"},
labels={
"background": 0,
"kidney": 1,
"tumor": 2
},
# regions_class_order=(1, 3, 2),
num_training_cases=210, file_ending='.nii.gz',
dataset_name=task_name, reference='none',
release='prerelease',
description="KiTS2019")
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('input_folder', type=str,
help="The downloaded and extracted KiTS2023 dataset (must have case_XXXXX subfolders)")
parser.add_argument('-d', required=False, type=int, default=209, help='nnU-Net Dataset ID, default: 220')
args = parser.parse_args()
amos_base = args.input_folder
convert_kits2023(amos_base, args.d)
# /root/autodl-tmp/kits
终端运行:
python Dataset210_KiTS2019.py /root/autodl-tmp/kits
总结: 这两种方法生成的文件夹不一样,所以dataset.json
文件也不一样。
3.4 数据预处理(开GPU)
3.4.1 重采样
在nnunet里面有重采样,但是冠军方法采用的是将所有病例的体素间距重采样为 3.22 x 1.62 x 1.62.
所以,可以自己进行采样(这里需要修改下nnunet里的重采样),然后再训练。nnUnet肾脏肿瘤分割实战(KiTS19)
我没有重采样,我用的是nnunet的预处理,我主要是看一个结果。这个代码我跑过,用CPU的时候运行时出现错误Killed,与下面的预处理遇到相同问题。所以我想应该开GPU跑就不会有问题,但是我没试过。
我添加了个reshaping的文件夹,是为了防止原文件出错。
import numpy as np
import SimpleITK as sitk
import os
'''
算法功能:进行重采样,将所有病例的体素间距重采样为 3.22 x 1.62 x 1.62.
代码出现错误:Killed.
'''
# 定义插值函数
def transform(image,newSpacing, resamplemethod=sitk.sitkNearestNeighbor):
# 设置一个Filter
resample = sitk.ResampleImageFilter()
# 初始的体素块尺寸
originSize = image.GetSize()
# 初始的体素间距
originSpacing = image.GetSpacing()
newSize = [
int(np.round(originSize[0] * originSpacing[0] / newSpacing[0])),
int(np.round(originSize[1] * originSpacing[1] / newSpacing[1])),
int(np.round(originSize[2] * originSpacing[2] / newSpacing[2]))
]
print('current size:',newSize)
# 沿着x,y,z,的spacing(3)
# The sampling grid of the output space is specified with the spacing along each dimension and the origin.
resample.SetOutputSpacing(newSpacing)
# 设置original
resample.SetOutputOrigin(image.GetOrigin())
# 设置方向
resample.SetOutputDirection(image.GetDirection())
resample.SetSize(newSize)
# 设置插值方式
resample.SetInterpolator(resamplemethod)
# 设置transform
resample.SetTransform(sitk.Euler3DTransform())
# 默认像素值 resample.SetDefaultPixelValue(image.GetPixelIDValue())
return resample.Execute(image)
# 给image进行插值,采用 B样条 插值
data_path = "/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_raw/Dataset209_KiTS2019/imagesTr"
# data_path311 = "/root/autodl-tmp/nnUNet/dataset/nnUNet_raw/Dataset040_KiTS/imagesTr"
data_path311 = "/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_raw/Dataset311_KiTS209/imagesTr"
for path in sorted(os.listdir(data_path)):
print(path)
img_path = os.path.join(data_path,path)
img_itk = sitk.ReadImage(img_path)
print('origin size:', img_itk.GetSize())
# image采用sitk.sitkBSpline插值
new_itk = transform(img_itk, [3.22, 1.62, 1.62], sitk.sitkBSpline) # sitk.sitkLinear
# sitk.WriteImage(new_itk, img_path)
data_path3 = os.path.join(data_path311,path)
sitk.WriteImage(new_itk, data_path3)
print('images is resampled!')
print('-'*20)
# 给mask进行插值,采用 NearestNeighbor 插值
label_path = "/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_raw/Dataset209_KiTS2019/labelsTr"
label_path311 = "/root/autodl-fs/nnUNet-wh/DATASET/nnUNet_raw/Dataset311_KiTS209/labelsTr"
for path in sorted(os.listdir(label_path)):
print(path)
img_path = os.path.join(label_path,path)
img_itk = sitk.ReadImage(img_path)
print('origin size:', img_itk.GetSize())
# segment采用sitk.sitkNearestNeighbor插值
new_itk = transform(img_itk, [3.22, 1.62, 1.62])
# sitk.WriteImage(new_itk, img_path)
label_path3 = os.path.join(label_path311,path)
sitk.WriteImage(new_itk, label_path3)
print('labels is resampled!')
# 测试代码哪里有问题,检查过了,还是出现了Killed的问题,所以我猜测应该是得开GPU。
# try:
# # Your existing code
# except Exception as e:
# print(f"An error occurred: {e}")
3.4.2 nnunet预处理
nnUNetv2_plan_and_preprocess -d 209 --verify_dataset_integrity
用CPU时,数据预处理出现错误,错误显示:Killed.
用GPU时,数据预处理正常显示:
四、训练
修改了epoch=500,我跑过一轮发现400左右就差不多了
autodl-fs/nnUNet-wh/nnunetv2/training/nnUNetTrainer/nnUNetTrainer.py
4.1 终端运行代码
209为数据文件夹编号,2d为模型(2d, 3d_fullres, 3d_lowres, 3d_cascade_fullres
) ,0
为五折交叉验证中的第0折(0-4),即210个数据分5份,其中168个数据用来训练,42个数据用来验证。all
是210个数据均用来训练,得到一个模型,没有验证。5
不是5折交叉验证,nnUNet会以4:1的比例随机选择训练集和验证集,来自nnUNetV2使用教程,超详细!!(使用MSD十项全能数据集)
4.1.1 一折一折的跑
nnUNetv2_train DATASET_NAME_OR_ID UNET_CONFIGURATION FOLD [additional options, see -h]
依次执行:
nnUNetv2_train 209 2d 0
# nnUNetv2_train 209 2d 1
# nnUNetv2_train 209 2d 2
# nnUNetv2_train 209 2d 3
# nnUNetv2_train 209 2d 4
# nnUNetv2_train 209 2d all
注: 如果您希望使用单个模型进行预测,请将all
fold进行训练,并在nnUNetv2_predict -f all
中使用的。
- 终端运行保存日志(我没试过)
nnUNetv2_train 209 2d 0 > train.log 2>&1
- 运行中断,继续运行
nnUNetv2_train 209 2d 0 --c
注: 训练后的结构在autodl-fs/nnUNet-wh/DATASET/nnUNet_trained_models/Dataset209_KiTS2019/nnUNetTrainer__nnUNetPlans__2d/fold_1
第0
折的验证集42个数据预测结果在 autodl-fs/nnUNet-wh/DATASET/nnUNet_trained_models/Dataset209_KiTS2019/nnUNetTrainer__nnUNetPlans__2d/fold_0/validation
4.1.2 自动运行5折
- 可以创建一个tst.sh 文件,并写入此代码:
for fold in {0..4}
do
# echo "nnUNetv2_train 1 3d_lowres $fold"
nnUNetv2_train 1 3d_lowres $fold
done
- 在虚拟环境下的终端运行此 tst.sh 文件,使用 source命令可以执行脚本(为什么在shell脚本中无法使用cd?问题原因及解决方法、nnU-Net v2的环境配置到训练自己的数据集),具体命令行为:
source /root/autodl-tmp/nnU-Net/sts.sh
4.2 如果预测时,想用五折验证的话,终端运行代码
这会让 nnU-Net 在最终验证期间
保存 softmax 输出。它们是必需的。导出的 softmax 预测非常大,因此可能占用大量磁盘空间,因此默认情况下不启用此功能。如果您最初没有使用--npz
标记运行,但现在需要 softmax 预测,请使用以下命令重新运行验证:
nnUNetv2_train DATASET_NAME_OR_ID UNET_CONFIGURATION FOLD --val --npz
依次执行:
nnUNetv2_train 209 2d 0 --val --npz
nnUNetv2_train 209 2d 1 --val --npz
nnUNetv2_train 209 2d 1 --val --npz
nnUNetv2_train 209 2d 1 --val --npz
nnUNetv2_train 209 2d 1 --val --npz
此处实际是把验证集42个数据重新跑下,保存 softmax 输出。
4.3 遇到的问题
运行一段时间,中断后,重新运行时,出现错误: ValueError: mmap length is greater than file size
解决方法:进入指定文件夹中autodl-fs/nnUNet-wh/DATASET/nnUNet_raw/Dataset209_KiTS2019/imagesTr
,执行(删除所有npy文件
)
rm *.npy
五、预测
5.1 单个模型测试
- 单个模型预测,下面例子
-f 0
是第0折模型。
nnUNetv2_predict -i INPUT_FOLDER -o OUTPUT_FOLDER -d DATASET_NAME_OR_ID -c CONFIGURATION --save_probabilities
nnUNetv2_predict -i ${nnUNet_raw}/Dataset209_KiTS2019/imagesTs -o output0 -d 209 -c 2d -f 0
# nnUNetv2_predict -i ${nnUNet_raw}/Dataset209_KiTS2019/imagesTs -o output1 -d 209 -c 2d -f 1
# nnUNetv2_predict -i ${nnUNet_raw}/Dataset209_KiTS2019/imagesTs -o output2 -d 209 -c 2d -f 2
# nnUNetv2_predict -i ${nnUNet_raw}/Dataset209_KiTS2019/imagesTs -o output3 -d 209 -c 2d -f 3
# nnUNetv2_predict -i ${nnUNet_raw}/Dataset209_KiTS2019/imagesTs -o output4 -d 209 -c 2d -f 4
# nnUNetv2_predict -i ${nnUNet_raw}/Dataset209_KiTS2019/imagesTs -o output_all -d 209 -c 2d -f all
5.2 寻找最佳配置
- 自动确定最佳配置
一旦所需的配置经过训练(完全5折交叉验证),您就可以告诉 nnU-Net 自动识别最适合您的组合:
nnU-Netv2在服务器上使用全流程(小白边踩坑边学习的记录)
nnUNetv2_find_best_configuration DATASET_NAME_OR_ID -c CONFIGURATIONS
nnUNetv2_find_best_configuration -h
查看参数
nnUNetv2_find_best_configuration 101 -f 0 1 2 3 4 -c 2d
5.3 在5.1和5.2的基础上,集成模型预测
- 集成模型,可以包括任意数量的文件夹,由
nnUNetv2_predict
生成的带npz的预测文件夹。
nnUNetv2_ensemble -i FOLDER1 FOLDER2 ... -o OUTPUT_FOLDER -np NUM_PROCESSES
5.4 后处理
nnUNetv2_apply_postprocessing -i FOLDER_WITH_PREDICTIONS -o OUTPUT_FOLDER --pp_pkl_file POSTPROCESSING_FILE -plans_json PLANS_FILE -dataset_json DATASET_JSON_FILE
5.5 结果压缩,得到predictions.zip压缩包
zip predictions.zip prediction_*.nii.gz
六、评价
6.1 nnunet自带评价
修改文件夹,运行,查看summary.json报告。
6.2 自己计算
二分类,label只有0和1
,三维nii
数据(如果是二维数据,需要给一下数据读取方式。)实现所有指标,并将结果保存为Excel
。【理论+实践】史上最全-论文中常用的图像分割评价指标-附完整代码
# 计算三维下各种指标
from __future__ import absolute_import, print_function
import pandas as pd
import GeodisTK
import numpy as np
from scipy import ndimage
# pixel accuracy
def binary_pa(s, g):
"""
calculate the pixel accuracy of two N-d volumes.
s: the segmentation volume of numpy array
g: the ground truth volume of numpy array
"""
pa = ((s == g).sum()) / g.size
return pa
# Dice evaluation
def binary_dice(s, g):
"""
calculate the Dice score of two N-d volumes.
s: the segmentation volume of numpy array
g: the ground truth volume of numpy array
"""
assert (len(s.shape) == len(g.shape))
prod = np.multiply(s, g)
s0 = prod.sum()
dice = (2.0 * s0 + 1e-10) / (s.sum() + g.sum() + 1e-10)
return dice
# IOU evaluation
def binary_iou(s, g):
assert (len(s.shape) == len(g.shape))
# 两者相乘值为1的部分为交集
intersecion = np.multiply(s, g)
# 两者相加,值大于0的部分为交集
union = np.asarray(s + g > 0, np.float32)
iou = intersecion.sum() / (union.sum() + 1e-10)
return iou
# Hausdorff and ASSD evaluation
def get_edge_points(img):
"""
get edge points of a binary segmentation result
"""
dim = len(img.shape)
if (dim == 2):
strt = ndimage.generate_binary_structure(2, 1)
else:
strt = ndimage.generate_binary_structure(3, 1) # 三维结构元素,与中心点相距1个像素点的都是邻域
ero = ndimage.morphology.binary_erosion(img, strt)
edge = np.asarray(img, np.uint8) - np.asarray(ero, np.uint8)
return edge
def binary_hausdorff95(s, g, spacing=None):
"""
get the hausdorff distance between a binary segmentation and the ground truth
inputs:
s: a 3D or 2D binary image for segmentation
g: a 2D or 2D binary image for ground truth
spacing: a list for image spacing, length should be 3 or 2
"""
s_edge = get_edge_points(s)
g_edge = get_edge_points(g)
image_dim = len(s.shape)
assert (image_dim == len(g.shape))
if (spacing == None):
spacing = [1.0] * image_dim
else:
assert (image_dim == len(spacing))
img = np.zeros_like(s)
if (image_dim == 2):
s_dis = GeodisTK.geodesic2d_raster_scan(img, s_edge, 0.0, 2)
g_dis = GeodisTK.geodesic2d_raster_scan(img, g_edge, 0.0, 2)
elif (image_dim == 3):
s_dis = GeodisTK.geodesic3d_raster_scan(img, s_edge, spacing, 0.0, 2)
g_dis = GeodisTK.geodesic3d_raster_scan(img, g_edge, spacing, 0.0, 2)
dist_list1 = s_dis[g_edge > 0]
dist_list1 = sorted(dist_list1)
dist1 = dist_list1[int(len(dist_list1) * 0.95)]
dist_list2 = g_dis[s_edge > 0]
dist_list2 = sorted(dist_list2)
dist2 = dist_list2[int(len(dist_list2) * 0.95)]
return max(dist1, dist2)
# 平均表面距离
def binary_assd(s, g, spacing=None):
"""
get the average symetric surface distance between a binary segmentation and the ground truth
inputs:
s: a 3D or 2D binary image for segmentation
g: a 2D or 2D binary image for ground truth
spacing: a list for image spacing, length should be 3 or 2
"""
s_edge = get_edge_points(s)
g_edge = get_edge_points(g)
image_dim = len(s.shape)
assert (image_dim == len(g.shape))
if (spacing == None):
spacing = [1.0] * image_dim
else:
assert (image_dim == len(spacing))
img = np.zeros_like(s)
if (image_dim == 2):
s_dis = GeodisTK.geodesic2d_raster_scan(img, s_edge, 0.0, 2)
g_dis = GeodisTK.geodesic2d_raster_scan(img, g_edge, 0.0, 2)
elif (image_dim == 3):
s_dis = GeodisTK.geodesic3d_raster_scan(img, s_edge, spacing, 0.0, 2)
g_dis = GeodisTK.geodesic3d_raster_scan(img, g_edge, spacing, 0.0, 2)
ns = s_edge.sum()
ng = g_edge.sum()
s_dis_g_edge = s_dis * g_edge
g_dis_s_edge = g_dis * s_edge
assd = (s_dis_g_edge.sum() + g_dis_s_edge.sum()) / (ns + ng)
return assd
# relative volume error evaluation
def binary_relative_volume_error(s_volume, g_volume):
s_v = float(s_volume.sum())
g_v = float(g_volume.sum())
assert (g_v > 0)
rve = abs(s_v - g_v) / g_v
return rve
def compute_class_sens_spec(pred, label):
"""
Compute sensitivity and specificity for a particular example
for a given class for binary.
Args:
pred (np.array): binary arrary of predictions, shape is
(height, width, depth).
label (np.array): binary array of labels, shape is
(height, width, depth).
Returns:
sensitivity (float): precision for given class_num.
specificity (float): recall for given class_num
"""
tp = np.sum((pred == 1) & (label == 1))
tn = np.sum((pred == 0) & (label == 0))
fp = np.sum((pred == 1) & (label == 0))
fn = np.sum((pred == 0) & (label == 1))
sensitivity = tp / (tp + fn)
specificity = tn / (tn + fp)
return sensitivity, specificity
def get_evaluation_score(s_volume, g_volume, spacing, metric):
if (len(s_volume.shape) == 4):
assert (s_volume.shape[0] == 1 and g_volume.shape[0] == 1)
s_volume = np.reshape(s_volume, s_volume.shape[1:])
g_volume = np.reshape(g_volume, g_volume.shape[1:])
if (s_volume.shape[0] == 1):
s_volume = np.reshape(s_volume, s_volume.shape[1:])
g_volume = np.reshape(g_volume, g_volume.shape[1:])
metric_lower = metric.lower()
if (metric_lower == "dice"):
score = binary_dice(s_volume, g_volume)
elif (metric_lower == "iou"):
score = binary_iou(s_volume, g_volume)
elif (metric_lower == 'assd'):
score = binary_assd(s_volume, g_volume, spacing)
elif (metric_lower == "hausdorff95"):
score = binary_hausdorff95(s_volume, g_volume, spacing)
elif (metric_lower == "rve"):
score = binary_relative_volume_error(s_volume, g_volume)
elif (metric_lower == "volume"):
voxel_size = 1.0
for dim in range(len(spacing)):
voxel_size = voxel_size * spacing[dim]
score = g_volume.sum() * voxel_size
else:
raise ValueError("unsupported evaluation metric: {0:}".format(metric))
return score
if __name__ == '__main__':
import os
import nibabel as nib
seg_path = '你的分割结果文件夹'
gd_path = '你的label文件夹'
save_dir = 'excel 存放文件夹'
seg = sorted(os.listdir(seg_path))
dices = []
hds = []
rves = []
case_name = []
senss = []
specs = []
for name in seg:
if not name.startswith('.') and name.endswith('nii.gz'):
# 加载label and segmentation image
seg_ = nib.load(os.path.join(seg_path, name))
seg_arr = seg_.get_fdata().astype('float32')
gd_ = nib.load(os.path.join(gd_path, name))
gd_arr = gd_.get_fdata().astype('float32')
case_name.append(name)
# 求hausdorff95距离
hd_score = get_evaluation_score(seg_arr, gd_arr, spacing=None, metric='hausdorff95')
hds.append(hd_score)
# 求体积相关误差
rve = get_evaluation_score(seg_arr, gd_arr, spacing=None, metric='rve')
rves.append(rve)
# 求dice
dice = get_evaluation_score(seg_.get_fdata(), gd_.get_fdata(), spacing=None, metric='dice')
dices.append(dice)
# 敏感度,特异性
sens, spec = compute_class_sens_spec(seg_.get_fdata(), gd_.get_fdata())
senss.append(sens)
specs.append(spec)
# 存入pandas
data = {'dice': dices, 'RVE': rves, 'Sens': senss, 'Spec': specs, 'HD95': hds}
df = pd.DataFrame(data=data, columns=['dice', 'RVE', 'Sens', 'Spec', 'HD95'], index=case_name)
df.to_csv(os.path.join(save_dir, 'metrics.csv'))
七、数据处理-excel快速实现"平均值±标准差"
- 选中插入公式的地方
- 插入函数,复制公式
=ROUND(AVERAGE(B2:B39),4) *100&"±"&ROUND(STDEV(B2:Bx), 4) *100
- 依次选中红框的数字,然后点number1,选中一列数据
- 得到计算的平均值和标准差,然后右拉,直接计算其他两列的值。