Technical dry goods | Hausdorff distance for image segmentation based on MindSpore

今天带来的内容是Hausdorff distance An introduction to the principle of Hausdorff distanceMindSpore的实现代码.

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

· Mean surface distance 平均表面距离

· Hausdorff distance 豪斯多夫距离（也被称为max_surface_distance 最大表面距离MSD）

· Surface overlap 表面重叠度

· Surface dice 表面dice值

· Volumetric dice 三维dice值

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

Hausdorff distance介绍

关于这个Hausdorff Distance Calculation of Hausdorff distance,There is really a lot of information online,But it feels a little messy,Almost all plagiarized from each other,讲的也不是很清楚,What's more, what many people have introduced is wrong.Here I will explain it according to my ideas（Of course, I am also the leader of a hundred families,进行一个总结,好吧,其实也是借鉴）.

之前的文章已经讲过Dicecoefficient（基于MindSporeComplete image segmentation efficiently,实现Dice！）,Dice对mask的内部填充比较敏感,而Hausdorff distance 对分割出的边界比较敏感,所以主要是用Hausdorff distanceUsed in image segmentation tasks.

Hausdorff distance原理及公式

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

这里,

· 式(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)两者中的较大者,它度量了两个点集间的最大不匹配程度.

Diagram and summary of calculation process

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

· 取A集合中的一点a0,计算a0到B集合中所有点的距离,保留最短的距离d0

· 遍历A集合中所有点,图中一共两点a0和a1,计算出d0和d1

· 比较所有的距离{ d0, d1 },选出最长的距离d1

· 这个最长的距离就是h,它是A→B的单向豪斯多夫距离,记为h( A, B )

· 对于A集合中任意一点a,我们可以确定,以点a为圆心,h为半径的圆内部必有B集合中的

· 交换A集合和B集合的角色,计算B→A的单向豪斯多夫距离h( B, A ),选出h( A, B )和h( B, A )中最长的距离,就是A,B集合的双向豪斯多夫距离

看完公式,Calculation process with diagrams and summaries,Are you asking for a lot,哈哈哈哈,Anyway, that's how I understand it.

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):
    # Generate a region of interest（ROI）.Clipping and spatial data are supported.The center and size of the space should be provided,如果没有,则必须提供ROIstart and end coordinates of .
    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距离的参数,Euclidean is supported、"chessboard"、"taxicab"Three measurement methods.
        self.distance_metric = validator.check_string(distance_metric, string_list, "distance_metric")
        # Represents the maximum distance quantile,取值范围为0-100,It represents computational steps4中,The selected distance covers the percentage of distance
        self.percentile = percentile if percentile is None else validator.check_value_type("percentile",
        # Can be divided into directional and non-directionalHausdorff距离,Default is non-directionalHausdorff距离,指定percentile参数得到HausdorffPercentile of distance.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):
        # Return inclusion by shortening or repeating the inputdim值的元组.
        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):
        # Check if it is a tuple
        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.percentileRepresents the maximum distance quantile,取值范围为0-100,It represents computational steps4中,The selected distance covers the percentage of distance
        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):
        # Use the Euclidean method to find surface distances
        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):
        """
        Calculate directional or non-directionalHausdorff 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)
        # Hausdorff of computational orientation
        if self.directed:
            return hd
        # Compute the non-directional Hausdorff
        hd2 = self._calculate_percent_hausdorff_distance(self.y_edges, self.y_pred_edges)
        # If the calculation is directional,那直接返回hd,If it is computing non-directional,那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 参数说明

Represents the maximum distance quantile,取值范围为0-100,It represents computational steps4中,The selected distance covers the percentage of distance,Usually selected95%,Then in the calculation step4is not the maximum distance selected in ,而是将距离从大到小排列后,取排名为5%的距离.这么做的目的是为了排除一些离群点所造成的不合理的距离,保持整体数值的稳定性.所以Hausdorff distance也被称为Hausdorff distance-95%.