图像分割评估指标之Hausdorff distance 豪斯多夫距离

我又来更新啦，今天带来的内容是Hausdorff distance 豪斯多夫距离的原理介绍及MindSpore的实现代码。

欢迎关注知乎： 世界是我改变的

知乎上的原文链接

当我们评价图像分割的质量和模型表现时，经常会用到各类表面距离的计算。比如：

• Mean surface distance 平均表面距离
• Hausdorff distance 豪斯多夫距离
• Surface overlap 表面重叠度
• Surface dice 表面dice值
• Volumetric dice 三维dice值
• max_surface_distance 最大表面距离

等等等等，都可以称为表面距离系列了。今天简单的讲解一下Hausdorff distance 豪斯多夫距离。

一. Hausdorff distance介绍

关于这个Hausdorff Distance 豪斯多夫距离的计算，网上资料真的好多好多，但感觉有些乱糟糟的，几乎都是互相抄袭，讲的也不是很清楚，更有甚者很多人介绍的是错的。这里我按照我的思路来讲解一下。（当然我这也是采百家之长，进行一个总结，好吧，其实也是借鉴。）

之前的文章已经讲过Dice系数了，Dice对mask的内部填充比较敏感，而Hausdorff distance 对分割出的边界比较敏感，所以主要是用Hausdorff distance使用在图像分割任务中。

1. Hausdorff distance原理及公式

Hausdorff distance是描述两组点集之间相似程度的一种量度，它是两个点集之间距离的一种定义形式：假设有两组集合A={a1,…,ap},B={b1,…,bq},则这两个点集合之间的Hausdorff distance定义为：

其中,

‖·‖是点集A和B点集间的距离范式。(如:L2或Euclidean距离).

这里,

• 式(1)称为双向Hausdorff distance，是Hausdorff distance的最基本形式；
• 式(2)中的h(A，B)和h(B，A)分别称为从A集合到B集合和从B集合到A集合的单向Hausdorff距离。即h(A,B)实际上首先对点集A中的每个点ai到距离此点ai最近的B集合中点bj之间的距离‖ai-bj‖进行排序，然后取该距离中的最大值作为h(A,B)的值；
• h(B,A)同理可得。
• 由式(1)知，双向Hausdorff距离H(A，B)是单向距离h(A，B)和h(B,A)两者中的较大者，它度量了两个点集间的最大不匹配程度。

2. 图示及计算流程总结

刚说了那么多，是不是也不是很清楚，只看公式确实是一件不好玩的事情，那我用网上常用的图来说明一下，还有一个比较简单清晰的计算流程。给定两个点集合A{ a0, a1, … }和B{ b0, b1, b2, …}

1. 取A集合中的一点a0，计算a0到B集合中所有点的距离，保留最短的距离d0
2. 遍历A集合中所有点，图中一共两点a0和a1，计算出d0和d1
3. 比较所有的距离{ d0, d1 }，选出最长的距离d1
4. 这个最长的距离就是h，它是A→B的单向豪斯多夫距离，记为h( A, B
5. 对于A集合中任意一点a，我们可以确定，以点a为圆心，h为半径的圆内部必有B集合中的
6. 交换A集合和B集合的角色，计算B→A的单向豪斯多夫距离h( B, A )，选出h( A, B )和h( B, A )中最长的距离，就是A,B集合的双向豪斯多夫距离

看完公式，配合着图和总结的计算流程，是不是就请出许多了，哈哈哈哈，反正我是这样看懂的。。。。

二. Hausdorff distance的MindSpore代码实现

好了，原理已经讲完，话不多说，我们开始上代码。使用的是MindSpore框架实现的代码。

• MindSpore代码实现

"""HausdorffDistance."""

from collections import abc
from abc import ABCMeta
from scipy.ndimage import morphology
import numpy as np
from mindspore.common.tensor import Tensor
from mindspore._checkparam import Validator as validator
from .metric import Metric


class _ROISpatialData(metaclass=ABCMeta):
    # 产生感兴趣区域（ROI）。支持裁剪和空间数据。应提供空间的中心和大小，如果没有，则必须提供ROI的起点和终点坐标。
    def __init__(self, roi_center=None, roi_size=None, roi_start=None, roi_end=None):

        if roi_center is not None and roi_size is not None:
            roi_center = np.asarray(roi_center, dtype=np.int16)
            roi_size = np.asarray(roi_size, dtype=np.int16)
            self.roi_start = np.maximum(roi_center - np.floor_divide(roi_size, 2), 0)
            self.roi_end = np.maximum(self.roi_start + roi_size, self.roi_start)
        else:
            if roi_start is None or roi_end is None:
                raise ValueError("Please provide the center coordinates, size or start coordinates and end coordinates"
                                 " of ROI.")
            # ROI起始坐标
            self.roi_start = np.maximum(np.asarray(roi_start, dtype=np.int16), 0)
            # ROI终点坐标
            self.roi_end = np.maximum(np.asarray(roi_end, dtype=np.int16), self.roi_start)

    def __call__(self, data):
        sd = min(len(self.roi_start), len(self.roi_end), len(data.shape[1:]))
        slices = [slice(None)] + [slice(s, e) for s, e in zip(self.roi_start[:sd], self.roi_end[:sd])]
        return data[tuple(slices)]


class HausdorffDistance(Metric):
    def __init__(self, distance_metric="euclidean", percentile=None, directed=False, crop=True):
        super(HausdorffDistance, self).__init__()
        string_list = ["euclidean", "chessboard", "taxicab"]
        distance_metric = validator.check_value_type("distance_metric", distance_metric, [str])
        # 计算Hausdorff距离的参数，支持欧氏、"chessboard"、"taxicab"三种测量方法。
        self.distance_metric = validator.check_string(distance_metric, string_list, "distance_metric")
        # 表示最大距离分位数，取值范围为0-100，它表示的是计算步骤4中，选取的距离能覆盖距离的百分比
        self.percentile = percentile if percentile is None else validator.check_value_type("percentile",
        # 可分为定向和非定向Hausdorff距离，默认为非定向Hausdorff距离，指定percentile参数得到Hausdorff距离的百分位数。                                                                                 percentile, [float])
        self.directed = directed if directed is None else validator.check_value_type("directed", directed, [bool])
        self.crop = crop if crop is None else validator.check_value_type("crop", crop, [bool])
        self.clear()

    def _is_tuple_rep(self, tup, dim):
        # 通过缩短或重复输入返回包含dim值的元组。
        result = None
        if not self._is_iterable_sequence(tup):
            result = (tup,) * dim
        elif len(tup) == dim:
            result = tuple(tup)

        if result is None:
            raise ValueError(f"Sequence must have length {dim}, but got {len(tup)}.")

        return result

    def _is_tuple(self, inputs):
        # 判断是否是一个元组
        if not self._is_iterable_sequence(inputs):
            inputs = (inputs,)

        return tuple(inputs)

    def _is_iterable_sequence(self, inputs):
        if isinstance(inputs, Tensor):
            return int(inputs.dim()) > 0
        return isinstance(inputs, abc.Iterable) and not isinstance(inputs, str)

    def _create_space_bounding_box(self, image, func=lambda x: x > 0, channel_indices=None, margin=0):
        data = image[[*(self._is_tuple(channel_indices))]] if channel_indices is not None else image
        data = np.any(func(data), axis=0)
        nonzero_idx = np.nonzero(data)
        margin = self._is_tuple_rep(margin, data.ndim)

        box_start = list()
        box_end = list()
        for i in range(data.ndim):
            if nonzero_idx[i].size <= 0:
                raise ValueError("did not find nonzero index at the spatial dim {}".format(i))
            box_start.append(max(0, np.min(nonzero_idx[i]) - margin[i]))
            box_end.append(min(data.shape[i], np.max(nonzero_idx[i]) + margin[i] + 1))
        return box_start, box_end

    def _calculate_percent_hausdorff_distance(self, y_pred_edges, y_edges):
        
        surface_distance = self._get_surface_distance(y_pred_edges, y_edges)

        if surface_distance.shape == (0,):
            return np.inf

        if not self.percentile:
            return surface_distance.max()
        # self.percentile表示最大距离分位数，取值范围为0-100，它表示的是计算步骤4中，选取的距离能覆盖距离的百分比
        if 0 <= self.percentile <= 100:
            return np.percentile(surface_distance, self.percentile)

        raise ValueError(f"percentile should be a value between 0 and 100, get {self.percentile}.")

    def _get_surface_distance(self, y_pred_edges, y_edges):
        # 使用欧式方法求表面距离
        if not np.any(y_pred_edges):
            return np.array([])

        if not np.any(y_edges):
            dis = np.inf * np.ones_like(y_edges)
        else:
            if self.distance_metric == "euclidean":
                dis = morphology.distance_transform_edt(~y_edges)
            elif self.distance_metric == "chessboard" or self.distance_metric == "taxicab":
                dis = morphology.distance_transform_cdt(~y_edges, metric=self.distance_metric)

        surface_distance = dis[y_pred_edges]

        return surface_distance

    def clear(self):
        """清楚历史数据"""
        self.y_pred_edges = 0
        self.y_edges = 0
        self._is_update = False

    def update(self, *inputs):
        """
        更新输入数据
        """
        if len(inputs) != 3:
            raise ValueError('HausdorffDistance need 3 inputs (y_pred, y, label), but got {}'.format(len(inputs)))
        y_pred = self._convert_data(inputs[0])
        y = self._convert_data(inputs[1])
        label_idx = inputs[2]

        if y_pred.size == 0 or y_pred.shape != y.shape:
            raise ValueError("Labelfields should have the same shape, but got {}, {}".format(y_pred.shape, y.shape))

        y_pred = (y_pred == label_idx) if y_pred.dtype is not bool else y_pred
        y = (y == label_idx) if y.dtype is not bool else y

        res1, res2 = None, None
        if self.crop:
            if not np.any(y_pred | y):
                res1 = np.zeros_like(y_pred)
                res2 = np.zeros_like(y)

            y_pred, y = np.expand_dims(y_pred, 0), np.expand_dims(y, 0)
            box_start, box_end = self._create_space_bounding_box(y_pred | y)
            cropper = _ROISpatialData(roi_start=box_start, roi_end=box_end)
            y_pred, y = np.squeeze(cropper(y_pred)), np.squeeze(cropper(y))

        self.y_pred_edges = morphology.binary_erosion(y_pred) ^ y_pred if res1 is None else res1
        self.y_edges = morphology.binary_erosion(y) ^ y if res2 is None else res2
        self._is_update = True

    def eval(self):
        """
        计算定向或者非定向的Hausdorff distance.
        """
        # 要先执行clear操作
        if self._is_update is False:
            raise RuntimeError('Call the update method before calling eval.')

        # 计算A到B的距离
        hd = self._calculate_percent_hausdorff_distance(self.y_pred_edges, self.y_edges)
        # 计算定向的豪斯多夫
        if self.directed:
            return hd
        # 计算非定向的豪斯多夫
        hd2 = self._calculate_percent_hausdorff_distance(self.y_edges, self.y_pred_edges)
        # 如果计算的是定向的，那直接返回hd，如果是计算非定向，那hd和hd2谁大返回谁
        return max(hd, hd2)
           

使用方法如下：

import numpy as np
from mindspore import Tensor
from mindspore.nn.metrics import HausdorffDistance

x = Tensor(np.array([[3, 0, 1], [1, 3, 0], [1, 0, 2]]))
y = Tensor(np.array([[0, 2, 1], [1, 2, 1], [0, 0, 1]]))
metric = HausdorffDistance()
metric.clear()
metric.update(x, y, 0)
distance = metric.eval()
print(distance)

1.4142135623730951
           

每个batch（比如两组数据）进行计算的时候如下：

import numpy as np
from mindspore import Tensor
from mindspore.nn.metrics import HausdorffDistance

x = Tensor(np.array([[3, 0, 1], [1, 3, 0], [1, 0, 2]]))
y = Tensor(np.array([[0, 2, 1], [1, 2, 1], [0, 0, 1]]))
metric = HausdorffDistance()
metric.clear()
metric.update(x, y, 0)

x1 = Tensor(np.array([[3, 0, 1], [1, 3, 0], [1, 0, 2]]))
y1 = Tensor(np.array([[0, 2, 1], [1, 2, 1], [0, 0, 1]]))
metric.update(x1, y1, 0)

distance = metric.eval()
print(distance)
           

• self.percentile参数说明

表示最大距离分位数，取值范围为0-100，它表示的是计算步骤4中，选取的距离能覆盖距离的百分比，一般是选取了95%，那么在计算步骤4中选取的不是最大距离，而是将距离从大到小排列后，取排名为5%的距离。这么做的目的是为了排除一些离群点所造成的不合理的距离，保持整体数值的稳定性。所以Hausdorff distance也被称为Hausdorff distance-95%。