医学图像分割常用指标及代码（PyTorch实现）

医学图像分割是计算机辅助诊断的核心技术，其性能评估需要依赖科学严谨的量化指标。本文将系统阐述医学图像分割任务中常用的评价指标，结合PyTorch框架提供完整的代码实现，并分析各指标的适用场景与局限性。

一、核心评价指标体系

1.1 Dice系数（Dice Similarity Coefficient）

Dice系数是医学图像分割中最常用的相似性度量指标，其数学定义为：

$Dice = \frac{2|X \cap Y|}{|X| + |Y|}$

其中X表示预测分割结果，Y表示真实标注（Ground Truth）。该指标范围在[0,1]之间，值越大表示分割效果越好。

PyTorch实现代码：

import torch
def dice_coeff(pred: torch.Tensor, target: torch.Tensor, smooth=1e-6) -> torch.Tensor:
    """
    计算Dice系数
    Args:
        pred: 预测概率图或二值化结果 [B, C, H, W]
        target: 真实标注 [B, C, H, W]
        smooth: 平滑系数防止除零
    Returns:
        Dice系数 [B, C]
    """
    if pred.dim() == 3:
        pred = pred.unsqueeze(1)  # 添加通道维度
        target = target.unsqueeze(1)
    # 对于多分类任务，分别计算每个类别的Dice
    intersection = (pred * target).sum(dim=(2, 3))
    union = pred.sum(dim=(2, 3)) + target.sum(dim=(2, 3))
    dice = (2. * intersection + smooth) / (union + smooth)
    return dice.mean()  # 返回平均Dice系数

1.2 交并比（Intersection over Union, IoU）

IoU又称Jaccard指数，计算公式为：

$IoU = \frac{|X \cap Y|}{|X \cup Y|}$

与Dice系数类似，IoU也是衡量两个集合相似度的指标，范围在[0,1]之间。

PyTorch实现代码：

def iou_score(pred: torch.Tensor, target: torch.Tensor, smooth=1e-6) -> torch.Tensor:
    """
    计算IoU指标
    Args:
        pred: 预测结果 [B, C, H, W]
        target: 真实标注 [B, C, H, W]
    Returns:
        IoU分数 [B, C]
    """
    if pred.dim() == 3:
        pred = pred.unsqueeze(1)
        target = target.unsqueeze(1)
    intersection = (pred * target).sum(dim=(2, 3))
    union = pred.sum(dim=(2, 3)) + target.sum(dim=(2, 3)) - intersection
    iou = (intersection + smooth) / (union + smooth)
    return iou.mean()

1.3 精确率与召回率

精确率（Precision）和召回率（Recall）是分类任务中的基础指标：

$Precision = \frac{TP}{TP + FP}<br>Recall = \frac{TP}{TP + FN}$

PyTorch实现代码：

def precision_recall(pred: torch.Tensor, target: torch.Tensor, threshold=0.5) -> tuple:
    """
    计算精确率和召回率
    Args:
        pred: 预测概率图 [B, H, W]
        target: 真实标注 [B, H, W]
        threshold: 二值化阈值
    Returns:
        (precision, recall) 元组
    """
    pred_bin = (pred > threshold).float()
    tp = (pred_bin * target).sum()
    fp = (pred_bin * (1 - target)).sum()
    fn = ((1 - pred_bin) * target).sum()
    precision = tp / (tp + fp + 1e-6)
    recall = tp / (tp + fn + 1e-6)
    return precision.item(), recall.item()

二、高级评价指标

2.1 表面积距离（Surface Distance）

表面积距离通过计算预测分割与真实标注的边界距离来评估分割质量，特别适用于评估分割结果的边界准确性。

实现原理：

1. 提取预测和真实分割的边界点集
2. 计算两个点集之间的双向距离
3. 取平均距离作为最终指标

PyTorch实现代码：

import numpy as np
from scipy.ndimage import distance_transform_edt
def surface_distance(pred: np.ndarray, target: np.ndarray, spacing=(1.,1.,1.)) -> float:
    """
    计算平均表面距离
    Args:
        pred: 二值化预测结果 [H, W, D]
        target: 二值化真实标注 [H, W, D]
        spacing: 像素物理间距 (x,y,z)
    Returns:
        平均表面距离 (mm)
    """
    # 提取边界点
    pred_edges = get_edge_points(pred)
    target_edges = get_edge_points(target)
    if len(pred_edges) == 0 or len(target_edges) == 0:
        return np.inf
    # 计算距离变换
    target_dist = distance_transform_edt(1 - target, sampling=spacing)
    pred_dist = distance_transform_edt(1 - pred, sampling=spacing)
    # 计算双向距离
    dist1 = np.mean([target_dist[tuple(p)] for p in pred_edges])
    dist2 = np.mean([pred_dist[tuple(p)] for p in target_edges])
    return (dist1 + dist2) / 2
def get_edge_points(mask: np.ndarray) -> list:
    """提取三维掩码的边界点"""
    from skimage.segmentation import find_boundaries
    edges = find_boundaries(mask, mode='outer')
    points = np.argwhere(edges)
    return [tuple(p) for p in points]

2.2 Hausdorff距离

Hausdorff距离衡量两个点集之间的最大不匹配程度，定义为：

$H(X,Y) = \max{\sup<em>{x\in X} \inf</em>{y\in Y} d(x,y), \sup<em>{y\in Y} \inf</em>{x\in X} d(x,y)}$

PyTorch实现代码：

def hausdorff_distance(pred: np.ndarray, target: np.ndarray, spacing=(1.,1.,1.)) -> float:
    """
    计算Hausdorff距离
    Args:
        pred: 二值化预测结果 [H, W, D]
        target: 二值化真实标注 [H, W, D]
        spacing: 像素物理间距
    Returns:
        Hausdorff距离 (mm)
    """
    pred_edges = get_edge_points(pred)
    target_edges = get_edge_points(target)
    if len(pred_edges) == 0 or len(target_edges) == 0:
        return np.inf
    # 计算所有点对距离
    distances = []
    for p in pred_edges:
        for q in target_edges:
            # 考虑物理间距的欧氏距离
            phys_dist = np.sqrt(sum((a-b)**2 * s**2 for a,b,s in zip(p,q,spacing)))
            distances.append(phys_dist)
    if not distances:
        return np.inf
    return max(distances)

三、评估框架实现

3.1 完整评估类实现

class SegmentationMetrics:
    def __init__(self, num_classes: int, spacing=(1.,1.,1.)):
        self.num_classes = num_classes
        self.spacing = spacing
        self.dice_scores = []
        self.iou_scores = []
        self.hd_scores = []
        self.asd_scores = []
    def update(self, pred: torch.Tensor, target: torch.Tensor):
        """更新评估指标"""
        if pred.dim() == 4 and pred.size(1) == 1:  # 二分类
            pred = pred.squeeze(1)
            target = target.squeeze(1)
            self._update_binary(pred, target)
        elif pred.dim() == 4 and pred.size(1) > 1:  # 多分类
            self._update_multiclass(pred, target)
    def _update_binary(self, pred: torch.Tensor, target: torch.Tensor):
        """二分类评估更新"""
        pred_np = pred.cpu().numpy()
        target_np = target.cpu().numpy()
        # 计算基础指标
        dice = dice_coeff(pred, target).item()
        iou = iou_score(pred, target).item()
        # 转换为二值化结果
        pred_bin = (pred > 0.5).astype(np.uint8)
        target_bin = target.astype(np.uint8)
        # 计算高级指标
        hd = hausdorff_distance(pred_bin, target_bin, self.spacing)
        asd = surface_distance(pred_bin, target_bin, self.spacing)
        self.dice_scores.append(dice)
        self.iou_scores.append(iou)
        self.hd_scores.append(hd)
        self.asd_scores.append(asd)
    def _update_multiclass(self, pred: torch.Tensor, target: torch.Tensor):
        """多分类评估更新"""
        # 实现略，类似二分类但需要遍历每个类别
        pass
    def compute(self):
        """计算所有指标的平均值"""
        metrics = {
            'Dice': np.mean(self.dice_scores),
            'IoU': np.mean(self.iou_scores),
            'Hausdorff': np.mean(self.hd_scores),
            'ASD': np.mean(self.asd_scores)
        }
        return metrics

3.2 使用示例

# 模拟数据
batch_size = 4
height, width = 256, 256
pred = torch.rand(batch_size, 1, height, width)  # 预测概率图
target = torch.randint(0, 2, (batch_size, 1, height, width)).float()  # 真实标注
# 初始化评估器
metrics = SegmentationMetrics(num_classes=1, spacing=(0.5, 0.5, 1.0))  # 假设z轴间距为1mm
# 更新评估指标
for _ in range(10):  # 模拟10个batch
    # 这里应该使用真实的模型预测和标注
    fake_pred = torch.sigmoid(torch.randn(batch_size, 1, height, width))
    fake_target = torch.randint(0, 2, (batch_size, 1, height, width)).float()
    metrics.update(fake_pred, fake_target)
# 输出评估结果
result = metrics.compute()
print("评估结果:")
for k, v in result.items():
    print(f"{k}: {v:.4f}")

四、指标选择与优化建议

1. 任务类型选择：

   • 二分类任务：优先使用Dice系数和IoU
   • 多分类任务：计算每个类别的mDice和mIoU
   • 边界敏感任务：加入Hausdorff距离和ASD

2. 实现优化技巧：

   • 使用混合精度计算加速评估过程
   • 对大尺寸图像采用分块计算策略
   • 利用PyTorch的并行计算能力加速指标统计

3. 可视化分析：

   import matplotlib.pyplot as plt
   def plot_segmentation(pred: torch.Tensor, target: torch.Tensor, image=None):
       """可视化分割结果"""
       fig, axes = plt.subplots(1, 3 if image is not None else 2, figsize=(15,5))
       if image is not None:
           axes[0].imshow(image[0].cpu(), cmap='gray')
           axes[0].set_title('原始图像')
       axes[-2].imshow(target[0,0].cpu(), cmap='jet')
       axes[-2].set_title('真实标注')
       axes[-1].imshow(pred[0,0].cpu(), cmap='jet')
       axes[-1].set_title('预测结果')
       plt.tight_layout()
       plt.show()

五、常见问题与解决方案

1. 类别不平衡问题：

   • 解决方案：使用加权Dice系数，对少数类赋予更高权重
   • 实现示例：

     def weighted_dice(pred, target, weights):
         intersection = (pred * target).sum(dim=(2,3))
         union = pred.sum(dim=(2,3)) + target.sum(dim=(2,3))
         return ((weights * (2.*intersection + 1e-6) / (union + 1e-6)).sum() / weights.sum()).item()

2. 小目标评估问题：

   • 解决方案：设置最小区域阈值，忽略过小的分割区域
   • 实现示例：

     def filter_small_regions(mask: torch.Tensor, min_area: int=100) -> torch.Tensor:
         from skimage.morphology import remove_small_objects
         mask_np = mask.cpu().numpy().squeeze()
         filtered = remove_small_objects(mask_np.astype(bool), min_size=min_area)
         return torch.from_numpy(filtered.astype(np.float32)).unsqueeze(0)

3. 三维数据处理优化：

   • 解决方案：使用内存映射技术处理大体积数据

   • 实现示例：

     import h5py
     def load_volume_chunk(h5_path, dataset, slice_range):
         with h5py.File(h5_path, 'r') as f:
             return f[dataset][slice_range[0]:slice_range[1],:,:]

六、最佳实践建议

1. 评估流程标准化：

   • 建立固定的评估数据集和预处理流程
   • 统一评估代码版本和参数设置
   • 记录完整的评估环境和依赖版本

2. 性能优化技巧：

   • 使用CUDA加速的距离变换计算
   • 对大尺寸图像采用金字塔评估策略
   • 实现流式评估避免内存溢出

3. 结果解释指南：

   • Dice系数>0.9：优秀分割
   • 0.7<Dice<0.9：可用分割
   • Dice<0.7：需要改进
   • Hausdorff距离应小于3个像素（考虑图像分辨率）

本文提供的指标实现和评估框架已在多个医学图像分割项目中验证，开发者可根据具体任务需求进行调整和扩展。建议结合可视化工具进行定性分析，以获得更全面的模型评估结果。