当前位置:   article > 正文

技术干货 | 基于 MindSpore 实现图像分割之豪斯多夫距离_图像分割surface distance

图像分割surface distance

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

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

· Mean surface distance 平均表面距离

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

· Surface overlap 表面重叠度

· Surface dice 表面dice值

· Volumetric dice 三维dice值

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

Hausdorff distance介绍

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

之前的文章已经讲过Dice系数了(基于MindSpore高效完成图像分割,实现Dice!),Dice对mask的内部填充比较敏感,而Hausdorff distance 对分割出的边界比较敏感,所以主要是用Hausdorff distance使用在图像分割任务中

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

图示及计算流程总结

刚说了那么多,是不是也不是很清楚,只看公式确实是一件不好玩的事情,那我用网上常用的图来说明一下,还有一个比较简单清晰的计算流程。给定两个点集合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集合的双向豪斯多夫距离

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

MindSpore 代码实现

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

  1. """HausdorffDistance."""
  2. from collections import abc
  3. from abc import ABCMeta
  4. from scipy.ndimage import morphology
  5. import numpy as np
  6. from mindspore.common.tensor import Tensor
  7. from mindspore._checkparam import Validator as validator
  8. from .metric import Metric
  9. class _ROISpatialData(metaclass=ABCMeta):
  10. # 产生感兴趣区域(ROI)。支持裁剪和空间数据。应提供空间的中心和大小,如果没有,则必须提供ROI的起点和终点坐标。
  11. def __init__(self, roi_center=None, roi_size=None, roi_start=None, roi_end=None):
  12. if roi_center is not None and roi_size is not None:
  13. roi_center = np.asarray(roi_center, dtype=np.int16)
  14. roi_size = np.asarray(roi_size, dtype=np.int16)
  15. self.roi_start = np.maximum(roi_center - np.floor_divide(roi_size, 2), 0)
  16. self.roi_end = np.maximum(self.roi_start + roi_size, self.roi_start)
  17. else:
  18. if roi_start is None or roi_end is None:
  19. raise ValueError("Please provide the center coordinates, size or start coordinates and end coordinates"
  20. " of ROI.")
  21. # ROI起始坐标
  22. self.roi_start = np.maximum(np.asarray(roi_start, dtype=np.int16), 0)
  23. # ROI终点坐标
  24. self.roi_end = np.maximum(np.asarray(roi_end, dtype=np.int16), self.roi_start)
  25. def __call__(self, data):
  26. sd = min(len(self.roi_start), len(self.roi_end), len(data.shape[1:]))
  27. slices = [slice(None)] + [slice(s, e) for s, e in zip(self.roi_start[:sd], self.roi_end[:sd])]
  28. return data[tuple(slices)]
  29. class HausdorffDistance(Metric):
  30. def __init__(self, distance_metric="euclidean", percentile=None, directed=False, crop=True):
  31. super(HausdorffDistance, self).__init__()
  32. string_list = ["euclidean", "chessboard", "taxicab"]
  33. distance_metric = validator.check_value_type("distance_metric", distance_metric, [str])
  34. # 计算Hausdorff距离的参数,支持欧氏、"chessboard"、"taxicab"三种测量方法。
  35. self.distance_metric = validator.check_string(distance_metric, string_list, "distance_metric")
  36. # 表示最大距离分位数,取值范围为0-100,它表示的是计算步骤4中,选取的距离能覆盖距离的百分比
  37. self.percentile = percentile if percentile is None else validator.check_value_type("percentile",
  38. # 可分为定向和非定向Hausdorff距离,默认为非定向Hausdorff距离,指定percentile参数得到Hausdorff距离的百分位数。percentile, [float])
  39. self.directed = directed if directed is None else validator.check_value_type("directed", directed, [bool])
  40. self.crop = crop if crop is None else validator.check_value_type("crop", crop, [bool])
  41. self.clear()
  42. def _is_tuple_rep(self, tup, dim):
  43. # 通过缩短或重复输入返回包含dim值的元组。
  44. result = None
  45. if not self._is_iterable_sequence(tup):
  46. result = (tup,) * dim
  47. elif len(tup) == dim:
  48. result = tuple(tup)
  49. if result is None:
  50. raise ValueError(f"Sequence must have length {dim}, but got {len(tup)}.")
  51. return result
  52. def _is_tuple(self, inputs):
  53. # 判断是否是一个元组
  54. if not self._is_iterable_sequence(inputs):
  55. inputs = (inputs,)
  56. return tuple(inputs)
  57. def _is_iterable_sequence(self, inputs):
  58. if isinstance(inputs, Tensor):
  59. return int(inputs.dim()) > 0
  60. return isinstance(inputs, abc.Iterable) and not isinstance(inputs, str)
  61. def _create_space_bounding_box(self, image, func=lambda x: x > 0, channel_indices=None, margin=0):
  62. data = image[[*(self._is_tuple(channel_indices))]] if channel_indices is not None else image
  63. data = np.any(func(data), axis=0)
  64. nonzero_idx = np.nonzero(data)
  65. margin = self._is_tuple_rep(margin, data.ndim)
  66. box_start = list()
  67. box_end = list()
  68. for i in range(data.ndim):
  69. if nonzero_idx[i].size <= 0:
  70. raise ValueError("did not find nonzero index at the spatial dim {}".format(i))
  71. box_start.append(max(0, np.min(nonzero_idx[i]) - margin[i]))
  72. box_end.append(min(data.shape[i], np.max(nonzero_idx[i]) + margin[i] + 1))
  73. return box_start, box_end
  74. def _calculate_percent_hausdorff_distance(self, y_pred_edges, y_edges):
  75. surface_distance = self._get_surface_distance(y_pred_edges, y_edges)
  76. if surface_distance.shape == (0,):
  77. return np.inf
  78. if not self.percentile:
  79. return surface_distance.max()
  80. # self.percentile表示最大距离分位数,取值范围为0-100,它表示的是计算步骤4中,选取的距离能覆盖距离的百分比
  81. if 0 <= self.percentile <= 100:
  82. return np.percentile(surface_distance, self.percentile)
  83. raise ValueError(f"percentile should be a value between 0 and 100, get {self.percentile}.")
  84. def _get_surface_distance(self, y_pred_edges, y_edges):
  85. # 使用欧式方法求表面距离
  86. if not np.any(y_pred_edges):
  87. return np.array([])
  88. if not np.any(y_edges):
  89. dis = np.inf * np.ones_like(y_edges)
  90. else:
  91. if self.distance_metric == "euclidean":
  92. dis = morphology.distance_transform_edt(~y_edges)
  93. elif self.distance_metric == "chessboard" or self.distance_metric == "taxicab":
  94. dis = morphology.distance_transform_cdt(~y_edges, metric=self.distance_metric)
  95. surface_distance = dis[y_pred_edges]
  96. return surface_distance
  97. def clear(self):
  98. """清除历史数据"""
  99. self.y_pred_edges = 0
  100. self.y_edges = 0
  101. self._is_update = False
  102. def update(self, *inputs):
  103. """
  104. 更新输入数据
  105. """
  106. if len(inputs) != 3:
  107. raise ValueError('HausdorffDistance need 3 inputs (y_pred, y, label), but got {}'.format(len(inputs)))
  108. y_pred = self._convert_data(inputs[0])
  109. y = self._convert_data(inputs[1])
  110. label_idx = inputs[2]
  111. if y_pred.size == 0 or y_pred.shape != y.shape:
  112. raise ValueError("Labelfields should have the same shape, but got {}, {}".format(y_pred.shape, y.shape))
  113. y_pred = (y_pred == label_idx) if y_pred.dtype is not bool else y_pred
  114. y = (y == label_idx) if y.dtype is not bool else y
  115. res1, res2 = None, None
  116. if self.crop:
  117. if not np.any(y_pred | y):
  118. res1 = np.zeros_like(y_pred)
  119. res2 = np.zeros_like(y)
  120. y_pred, y = np.expand_dims(y_pred, 0), np.expand_dims(y, 0)
  121. box_start, box_end = self._create_space_bounding_box(y_pred | y)
  122. cropper = _ROISpatialData(roi_start=box_start, roi_end=box_end)
  123. y_pred, y = np.squeeze(cropper(y_pred)), np.squeeze(cropper(y))
  124. self.y_pred_edges = morphology.binary_erosion(y_pred) ^ y_pred if res1 is None else res1
  125. self.y_edges = morphology.binary_erosion(y) ^ y if res2 is None else res2
  126. self._is_update = True
  127. def eval(self):
  128. """
  129. 计算定向或者非定向的Hausdorff distance.
  130. """
  131. # 要先执行clear操作
  132. if self._is_update is False:
  133. raise RuntimeError('Call the update method before calling eval.')
  134. # 计算A到B的距离
  135. hd = self._calculate_percent_hausdorff_distance(self.y_pred_edges, self.y_edges)
  136. # 计算定向的豪斯多夫
  137. if self.directed:
  138. return hd
  139. # 计算非定向的豪斯多夫
  140. hd2 = self._calculate_percent_hausdorff_distance(self.y_edges, self.y_pred_edges)
  141. # 如果计算的是定向的,那直接返回hd,如果是计算非定向,那hd和hd2谁大返回谁
  142. return max(hd, hd2)

使用方法如下:

  1. import numpy as np
  2. from mindspore import Tensor
  3. from mindspore.nn.metrics import HausdorffDistance
  4. x = Tensor(np.array([[3, 0, 1], [1, 3, 0], [1, 0, 2]]))
  5. y = Tensor(np.array([[0, 2, 1], [1, 2, 1], [0, 0, 1]]))
  6. metric = HausdorffDistance()
  7. metric.clear()
  8. metric.update(x, y, 0)
  9. distance = metric.eval()
  10. print(distance)
  11. 1.4142135623730951

每个batch(比如两组数据)进行计算的时候如下:

  1. import numpy as np
  2. from mindspore import Tensor
  3. from mindspore.nn.metrics import HausdorffDistance
  4. x = Tensor(np.array([[3, 0, 1], [1, 3, 0], [1, 0, 2]]))
  5. y = Tensor(np.array([[0, 2, 1], [1, 2, 1], [0, 0, 1]]))
  6. metric = HausdorffDistance()
  7. metric.clear()
  8. metric.update(x, y, 0)
  9. x1 = Tensor(np.array([[3, 0, 1], [1, 3, 0], [1, 0, 2]]))
  10. y1 = Tensor(np.array([[0, 2, 1], [1, 2, 1], [0, 0, 1]]))
  11. metric.update(x1, y1, 0)
  12. distance = metric.eval()
  13. print(distance)

self.percentile 参数说明

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

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/我家自动化/article/detail/793955
推荐阅读
  

闽ICP备14008679号