nnUnet肾脏肿瘤分割实战(KiTS19)
KiTS 数据集是一个用于医学影像分割的开源数据集，包含来自多个患者的CT扫描图像和对应的标签信息。该数据集经过预处理后存储在nnUNetrawdatabase/nnUNetrawdata/Task040KiTS目录下，包含训练集图像、训练集标签、验证集图像和验证集标签。预处理步骤包括裁剪（removing black borders）和重采样（resampling），以确保所有样本的空间分辨率一致。CT模态的体素间距被设置为3.22 x 1.62 x 1.62 mm³，并使用B样条插值进行插值操作以减少空间分辨率的变化。此外，还对标签中的错误进行了修正，并提供了完整的分割标签文件夹labelsTr。模型训练使用了nnUNet框架，并通过PyTorch实现 Dice 系数作为评价指标来评估分割性能。
"0": "CT",}
json_dict[‘labels’] = {
“0”: “background”,
“1”: “Kidney”,
“2”: “Tumor”
}
json_dict[‘numTraining’] = 210
json_dict[‘numTest’] = 90
json_dict[‘training’] = [{‘image’: “./imagesTr/%s.nii.gz” % i, “label”: “./labelsTr/%s.nii.gz” % i} for i in
cases[:210]]
json_dict[‘test’] = [“./imagesTs/%s.nii.gz” % i for i in
cases[210:]]
save_json(json_dict, os.path.join(out, “dataset.json”))
>
> 这里只是对数据集进行一个拷贝和重命名，不对原始数据进行修改。
>
>
>
运行代码后，整理好的数据集结构如下：该存储位置为nnUNet_raw_data_base下的子目录nnUNet_raw_data中的Task 49 KiTS。其中包含以下主要文件：
- 数据集信息记录在dataset.json中；
- 图像数据集位于imagesTr目录下；
- 包括多个病例的数据样本如case_... .nii.gz等。
├── imagesTs
│ ├── case_00210_0000.nii.gz
│ ├── case_00211_0000.nii.gz
│ ├── …
├── labelsTr
│ ├── case_00000.nii.gz
│ ├── case_00001.nii.gz
│ ├── …
>
> `dataset.json`文件保存了训练集图像、训练集标签、测试集图像等信息。
>
>
>
预处理阶段会根据`dataset.json`读取图像，如果想要剔除某个病例，直接在`dataset.json`修改就好。{
“description”: “kidney and kidney tumor segmentation”,
“labels”: {
“0”: “background”,
“1”: “Kidney”,
“2”: “Tumor”
},
“licence”: “”,
“modality”: {
“0”: “CT”
},
“name”: “KiTS”,
“numTest”: 90,
“numTraining”: 210,
“reference”: “KiTS data for nnunet”,
“release”: “0.0”,
“tensorImageSize”: “4D”,
“test”: [
“./imagesTs/case_00210.nii.gz”,
“./imagesTs/case_00211.nii.gz”,
…
],
“training”: [
{
“image”: “./imagesTr/case_00000.nii.gz”,
“label”: “./labelsTr/case_00000.nii.gz”
},
{
“image”: “./imagesTr/case_00001.nii.gz”,
“label”: “./labelsTr/case_00001.nii.gz”
},
…
]
}
提前准备三个文件夹，分别存放数据集、预处理数据和训练结果，配置好环境变量，具体细节可以参考我的[第一篇博文]()。
---
#### 2.数据预处理
nnUnet可以读取CT图像的模态信息、体素间距、灰度分布，自动进行重采样、裁剪以及归一化。
![在这里插入图片描述]()
nnUnet图像分割的自动方法配置(https://www.nature.com/articles/s41592-020-01008-z)
##### 重采样
不同时期，不同仪器的CT扫描仪，采样得到的CT图像具有不同的空间分辨率，重采样的目的是将所有的病例采样到相同的空间分辨率（体素间距）。
nnUnet的数据预处理`preprocess`自带重采样，但我试过两次之后效果并不好，重采样之后的图像尺寸太大了，于是我按照冠军论文里的方法自己写了个重采样，将所有病例的体素间距重采样为 `3.22 x 1.62 x 1.62`.
另外，论文中有提到case15和case37标签的错误，本来打算去掉，不过后来我去[KiTS19的github官网]()看了一下，官方已经作了修正。import numpy as np
import SimpleITK as sitk
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) 注意重采样的插值方法，我试过`SimpleITK`自带的多种插值方法，线性插值,三次插值以及B样条，比较发现B样条的效果是最好的。
因此，`image`采用`sitk.sitkBSpline`插值，`segment`采用`sitk.sitkNearestNeighbor`插值。
如果感兴趣可以自己尝试一下不同的插值方法，或者使用scipy等其他工具包进行重采样。data_path_ = "\\\\\\\\\\\\\\"/root/data/nnunet_raw_data_base/nnunet_raw_data
for name in sorted(os.listdir(data_dir)):
print(name)
image_dir = os.path.join(data_dir, name)
image_itk = sitk.ReadImage(image_dir)
print('origin size:', img_itk.GetSize())
new_itk = transform(img_itk, [3.22, 1.62, 1.62], sitk.sitkBSpline) # sitk.sitkLinear
sitk.WriteImage(new_itk, img_path)print(‘images is resampled!’)
print(‘-’*20)
labeled_path 包含子目录如：根目录下的数据集 nnUNet_base 中的 KiTS 任务的标签文件夹
遍历sorted(os.listdir(label_path))中的每个文件名path
输出当前处理的文件路径
构建完整的路径名
使用SimpleITK模块读取对应的图像数据
print('origin size:', img_itk.GetSize())
new_itk = transform(img_itk, [3.22, 1.62, 1.62])
sitk.WriteImage(new_itk, img_path)print(‘labels is resampled!’)
---
下面开始介绍nnUnet的数据预处理方法：
**输入指令：**python nnunet/experiment_planning/nnUNet_plan_and_preprocess.py -t 40 --verify_dataset_integrity
`verify_dataset_integrity`这里不再赘述，主要是根据验证数据集结构，第一次运行的时候最好还是加上。
##### 裁剪
裁剪的目的是裁去黑边，减少像素值为0的边缘区域，裁剪的时候保持空间分辨率等信息不变。此函数用于设置KiTS数据集的切片任务参数：task_string指定任务名称；override布尔值默认为False；num_threads设置多线程数量，默认采用系统默认值
该函数通过以下步骤操作：首先构建输出目录路径并将其赋值给变量cropped_out_dir；随后调用join函数将当前任务字符串与预设路径组合生成完整输出目录
为了确保创建父目录时文件夹已存在必要性而使用了maybe_mkdir_p辅助函数
if override and isdir(cropped_out_dir):
shutil.rmtree(cropped_out_dir)
maybe_mkdir_p(cropped_out_dir)
splitted_4d_output_dir_task = join(nnUNet_raw_data, task_string)
lists, _ = create_lists_from_splitted_dataset(splitted_4d_output_dir_task) # 创建裁剪列表
imgcrop = ImageCropper(num_threads, cropped_out_dir)
imgcrop.run_cropping(lists, overwrite_existing=override)
shutil.copy(join(nnUNet_raw_data, task_string, "dataset.json"), cropped_out_dir) **create\_lists\_from\_splitted\_dataset**加载所有的训练集的图像地址，lists一共有210个元素，每个元素包含图像和标签。def create_lists_from_splitted_dataset(base_folder_splitted):
lists = []
json_file = join(base_folder_splitted, "dataset.json")
with open(json_file) as jsn:
d = json.load(jsn)
training_files = d['training']
num_modalities = len(d['modality'].keys())
for tr in training_files:
cur_pat = []
for mod in range(num_modalities):
cur_pat.append(join(base_folder_splitted, "imagesTr", tr['image'].split("/")[-1][:-7] +
"\_%04.0d.nii.gz" % mod))
cur_pat.append(join(base_folder_splitted, "labelsTr", tr['label'].split("/")[-1]))
lists.append(cur_pat)
return lists, {int(i): d['modality'][str(i)] for i in d['modality'].keys()} 重点是这两个函数： imgcrop = ImageCropper(num_threads, cropped_out_dir)
imgcrop.run_cropping(lists, overwrite_existing=override) **ImageCropper**是一个类，包含10个方法。
重点是crop和run\_cropping两个方法:
* crop：裁剪到非零区域，返回data, seg, properties
* run\_cropping：执行裁剪操作，并且将结果保存为.npz文件(包含data和seg)，将size, spacing, origin, classes, size\_after\_cropping 等属性保存在.pkl文件。
![在这里插入图片描述]()
但是执行代码时，发现裁剪前后尺寸没有变化，可能是因为图像没有什么黑边 # 裁剪的时候seg!=None
def crop(data, properties, seg=None):
shape_before = data.shape # 原始尺寸
data, seg, bbox = crop_to_nonzero(data, seg, nonzero_label=-1) # 裁剪结果
shape_after = data.shape # 裁剪尺寸
print("before crop:", shape_before, "after crop:", shape_after, "spacing:",
np.array(properties["original\_spacing"]), "\n")
properties["crop\_bbox"] = bbox
properties['classes'] = np.unique(seg)
seg[seg < -1] = 0
properties["size\_after\_cropping"] = data[0].shape
return data, seg, properties ##### 数据分析
收集上一步裁剪得到的图像信息（尺寸、体素间距、灰度分布），为当前任务制定合适的训练计划(plan) # '/root/data/nnUNet\_raw\_data\_base/nnUNet\_cropped\_data/Task040\_KiTS'
cropped_out_dir = os.path.join(nnUNet_cropped_data, t)
# '/root/data/nnUNet\_preprocessed/Task040\_KiTS'
preprocessing_output_dir_this_task = os.path.join(preprocessing_output_dir, t)
# we need to figure out if we need the intensity propoerties. We collect them only if one of the modalities is CT
dataset_json = load_json(join(cropped_out_dir, 'dataset.json'))
modalities = list(dataset_json["modality"].values())
collect_intensityproperties = True if (("CT" in modalities) or ("ct" in modalities)) else False
dataset_analyzer = DatasetAnalyzer(cropped_out_dir, overwrite=False, num_processes=tf) # this class creates the fingerprint
_ = dataset_analyzer.analyze_dataset(collect_intensityproperties) # this will write output files that will be used by the ExperimentPlanner
maybe_mkdir_p(preprocessing_output_dir_this_task)
shutil.copy(join(cropped_out_dir, "dataset\_properties.pkl"), preprocessing_output_dir_this_task)
shutil.copy(join(nnUNet_raw_data, t, "dataset.json"), preprocessing_output_dir_this_task) 分析得到的`dataset_properties.pkl`结果如下：
![在这里插入图片描述]()
##### 创建数据指纹
根据上一步得到的数据集信息，针对不同的训练任务，制定合适的训练计划(plan) if planner_3d is not None:
if args.overwrite_plans is not None:
assert args.overwrite_plans_identifier is not None, "You need to specify -overwrite\_plans\_identifier"
exp_planner = planner_3d(cropped_out_dir, preprocessing_output_dir_this_task, args.overwrite_plans,
args.overwrite_plans_identifier)
else:
exp_planner = planner_3d(cropped_out_dir, preprocessing_output_dir_this_task)
exp_planner.plan_experiment()
if not dont_run_preprocessing: # double negative, yooo
exp_planner.run_preprocessing(threads)
if planner_2d is not None:
exp_planner = planner_2d(cropped_out_dir, preprocessing_output_dir_this_task)
exp_planner.plan_experiment()
if not dont_run_preprocessing: # double negative, yooo
exp_planner.run_preprocessing(threads) 预处理执行完毕，得到如下处理结果：位于nnUNet预处理文件夹内，
├── Task 40 KiTS数据集
├── annotation.json
├── annotation_info.pkl
├── groundtruth_segmentations
├── nnUNetData plans v2.1二维分割计划阶段 0
├── nnUNET plans v2-1二维分割计划阶段 0
└── validation_splits.pkl
这里生成的文件都可以打开来看看，对预处理方法和数据指纹有一个了解
* dataset.json在**数据获取**阶段产生
* daset\_properties为数据的 size, spacing, origin, classes, size\_after\_cropping 等属性
* gt\_segmentations为图像分割标签
* nnUNetData\_plans\_v2.1\_2D\_stage0和nnUNetData\_plans\_v2.1\_stage0是预处理后的数据集
* splits\_final.pkl是五折交叉验证划分的结果，一共210个病人，42为一折
* nnUNetPlansv2.1\_plans\*.pkl为训练计划，参考官方文档中的**edit\_plans\_files.md**可进行编辑
以`nnUNetPlansv2.1_plans_3D.pkl`为例，
![在这里插入图片描述]()
#### 3.模型训练
一行代码开始训练，执行过程以及调参可以参考我的博客[nnUnet代码解读–模型训练]()python nnunet/run/run_training.py CONFIGURATION TRAINER_CLASS_NAME TASK_NAME_OR_ID FOLD # FORMATTING
python nnunet/run/run_training.py 3d_fullres nnUNetTrainerV2 40 1
训练开始后，训练日志和训练结果记录在nnUNet\_trained\_models/nnUNet/3d\_fullres/Task040\_KiTS文件夹下UNetTrainer__nnUNetPlansv2.1
├── fold_1
│ ├── debug.json
│ ├── model_best.model
│ ├── model_best.model.pkl
│ ├── model_final_checkpoint.model
│ ├── model_final_checkpoint.model.pkl
│ ├── postprocessing.json
│ ├── progress.png
│ ├── training_log_2022_5_4_12_06_14.txt
│ ├── training_log_2022_5_5_10_30_05.txt
│ ├── validation_raw
│ └── validation_raw_postprocessed
![在这里插入图片描述]()
训练过程Loss曲线以及在线计算的Dice曲线
>
> 这里我想补充一下nnUnet的评价指标
>
>
>
**在线评价**
下面这段代码是nnUnet计算dice值的方法
先对每张图像中的每个类别分别计算tp, fp, fn，再对一个batch内的所有图像的tp, fp, fn求和，同时对一个batch求diceimport numpy as np
import torch
计算张量沿指定轴的总和。具体而言：
- 首先提取并转换输入轴为整数类型；
- 如果参数keepdim设为True，则会将输入沿着所有指定轴进行求和运算，并保留原始维度；
- 否则会按照降序排列指定轴并依次进行求和运算；
- 最终返回计算后的结果张量。
def run_online_evaluation(output, target):
torch.Size([b,num_classes, 80, 160, 160]) torch.Size([b,1, 80, 160, 160])
with torch.no_grad():
num_classes = output.shape[1]
output_softmax = torch.softmax(output,dim=1)
output_seg = output_softmax.argmax(1) # [b,80,160,160]
target = target[:, 0] # [b,80,160,160]
axes = tuple(range(1, len(target.shape))) # (1,2,…,num_classes)
tp_hard = torch.zeros((target.shape[0], num_classes - 1)).to(output_seg.device.index) # [b,num_classes-1]
fp_hard = torch.zeros((target.shape[0], num_classes - 1)).to(output_seg.device.index) # [b,num_classes-1]
fn_hard = torch.zeros((target.shape[0], num_classes - 1)).to(output_seg.device.index) # [b,num_classes-1]
for c in range(1, num_classes):
tp_hard[:, c - 1] = sum_tensor((output_seg == c).float() * (target == c).float(), axes=axes)
fp_hard[:, c - 1] = sum_tensor((output_seg == c).float() * (target != c).float(), axes=axes)
fn_hard[:, c - 1] = sum_tensor((output_seg != c).float() * (target == c).float(), axes=axes)
[b,num_classes-1] -> [num_classes-1,]
tp_hard = tp_hard.sum(0, keepdim=False).detach().cpu().numpy()
fp_hard = fp_hard.sum(0, keepdim=False).detach().cpu().numpy()
fn_hard = fn_hard.sum(0, keepdim=False).detach().cpu().numpy()
print(list((2 \* tp_hard) / (2 \* tp_hard + fp_hard + fn_hard + 1e-8)))
print(list(tp_hard))
print(list(fp_hard))
print(list(fn_hard))In the main module execution,
outputs are initialized as a random tensor with dimensions [batch_size=4,
channels=3,
height=80,
width=160].
Targets are generated as random integers within [low=0,
num_classes=3],
shaped to match the output dimensions.
The online evaluation process is then executed using these outputs and targets.
>
> 但是我觉得直接对一个batch累加求dice不够准确，因为不同图像的目标区域大小不同，目标区域大的图像对目标区域小的图像影响太大了。
>
>
>
>
> 比较好的评价方法是应该对batch内的每张图像分别求dice，然后求平均。
>
>
>
下面这段代码中，作者也提到：
训练过程中的在线评价，只是对Dice值的一个估计，并不能代表最终的dice.
整体思路就是把每个batch当做一张图像去求的dice，迭代一个epoch之后，再对每个batch的dice求平均。
验证时，每个epoch中batch的数量取决于`num_val_batches_per_epoch` def finish\_online\_evaluation(self):
self.online_eval_tp = np.sum(self.online_eval_tp, 0)
self.online_eval_fp = np.sum(self.online_eval_fp, 0)
self.online_eval_fn = np.sum(self.online_eval_fn, 0)
global_dc_per_class = [i for i in [2 \* i / (2 \* i + j + k) for i, j, k in
zip(self.online_eval_tp, self.online_eval_fp, self.online_eval_fn)]
if not np.isnan(i)]
self.all_val_eval_metrics.append(np.mean(global_dc_per_class))
self.print_to_log_file("Average global foreground Dice:", [np.round(i, 4) for i in global_dc_per_class])
网上资源浩如烟海，在这种情况下如果不能实现知识体系化管理，在遇到技术难题时往往停留在表面认识无法深入挖掘其内在规律就很难实现真正的技术进步。
单个个体的能力或许较强...然而集体的力量往往能够取得更大的成就...不论你是正从事IT行业的资深从业者还是抱着对IT行业的浓厚兴趣的新手入 join our community（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们共同进步与成长！
val_fn）
如果i不是NaN值，则执行以下操作：
将计算得到的全局分类性能指标平均值添加到all_val_eval_metrics列表中。
self.print_to_log_file("Average global foreground Dice:", [np.round(i, 4) for i in global_dc_per_class])[外链图片转存中…(img-BRzbGQ6E-1714240581584)]
[外链图片转存中…(img-IyANQtLc-1714240581585)]
网上充斥着海量的学习资料，在知识获取方面确实让人应接不暇。然而，在缺乏系统性的知识架构下，在面对技术难题时往往停留在表面探讨而非深入研究与彻底解决，在这样的情况下难以实现真正的技术进步与提升。
一个人可以走得效率高，但是一群人才能走得更远！无论你是正在从事IT行业的老鸟还是对IT行业充满热情的新手，我们都欢迎你加入我们的圈子（技术交流；学习资源；职场吐槽；大厂内推；面试辅导），让我们一起学习成长！