圖像分割評估名額之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%。