scikit-image · jni · Jul 2, 2020 · Sep 29, 2013 · Sep 29, 2013 · Sep 29, 2013
diff --git a/benchmarks/benchmark_metrics.py b/benchmarks/benchmark_metrics.py
@@ -0,0 +1,18 @@
+import numpy as np
+
+from skimage import metrics
+
+
+class SetMetricsSuite(object):
+    shape = (6, 6)
+    coords_a = np.zeros(shape, dtype=np.bool)
+    coords_b = np.zeros(shape, dtype=np.bool)
+
+    def setup(self):
+        points_a = (1, 0)
+        points_b = (5, 2)
+        self.coords_a[points_a] = True
+        self.coords_b[points_b] = True
+
+    def time_hausdorff(self):
+        metrics.hausdorff_distance(self.coords_a, self.coords_b)
diff --git a/doc/examples/segmentation/plot_hausdorff_distance.py b/doc/examples/segmentation/plot_hausdorff_distance.py
@@ -0,0 +1,74 @@
+"""
+==================
+Hausdorff Distance
+==================
+
+This example shows how to calculate the Hausdorff distance between two sets of
+points. The `Hausdorff distance
+<https://en.wikipedia.org/wiki/Hausdorff_distance>`__ is the maximum distance
+between any point on the first set and its nearest point on the second set,
+and vice-versa.
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+from skimage import metrics
+
+shape = (60, 60)
+image = np.zeros(shape)
+
+# Create a diamond-like shape where the four corners form the 1st set of points
+x_diamond = 30
+y_diamond = 30
+r = 10
+
+fig, ax = plt.subplots()
+plt_x = [0, 1, 0, -1]
+plt_y = [1, 0, -1, 0]
+
+set_ax = [(x_diamond + r * x) for x in plt_x]
+set_ay = [(y_diamond + r * y) for y in plt_y]
+plt.plot(set_ax, set_ay, 'or')
+
+# Create a kite-like shape where the four corners form the 2nd set of points
+x_kite = 30
+y_kite = 30
+x_r = 15
+y_r = 20
+
+set_bx = [(x_kite + x_r * x) for x in plt_x]
+set_by = [(y_kite + y_r * y) for y in plt_y]
+plt.plot(set_bx, set_by, 'og')
+
+# Set up the data to compute the hausdorff distance
+coords_a = np.zeros(shape, dtype=np.bool)
+coords_b = np.zeros(shape, dtype=np.bool)
+for x, y in zip(set_ax, set_ay):
+    coords_a[(x, y)] = True
+
+for x, y in zip(set_bx, set_by):
+    coords_b[(x, y)] = True
+
+# Call the hausdorff function on the coordinates
+metrics.hausdorff_distance(coords_a, coords_b)
+
+# Plot the lines that shows the length of the hausdorff distance
+x_line = [30, 30]
+y_line = [20, 10]
+plt.plot(x_line, y_line, 'y')
+
+x_line = [30, 30]
+y_line = [40, 50]
+plt.plot(x_line, y_line, 'y')
+
+# Plot circles to show that at this distance, the hausdorff distance can
+# travel to its nearest neighbor (in this case, from the kite to diamond)
+ax.add_artist(plt.Circle((30, 10), 10, color='y', fill=None))
+ax.add_artist(plt.Circle((30, 50), 10, color='y', fill=None))
+ax.add_artist(plt.Circle((15, 30), 10, color='y', fill=None))
+ax.add_artist(plt.Circle((45, 30), 10, color='y', fill=None))
+
+ax.imshow(image, cmap=plt.cm.gray)
+ax.axis((0, 60, 60, 0))
+plt.show()
diff --git a/skimage/metrics/__init__.py b/skimage/metrics/__init__.py
@@ -5,12 +5,14 @@
                              normalized_root_mse,
                              peak_signal_noise_ratio)
 from ._structural_similarity import structural_similarity
+from .set_metrics import hausdorff_distance
 
 __all__ = ['adapted_rand_error',
            'variation_of_information',
            'contingency_table',
            'mean_squared_error',
            'normalized_root_mse',
            'peak_signal_noise_ratio',
-           'structural_similarity'
+           'structural_similarity',
+           'hausdorff_distance'
            ]
diff --git a/skimage/metrics/set_metrics.py b/skimage/metrics/set_metrics.py
@@ -0,0 +1,51 @@
+import numpy as np
+from scipy.spatial import cKDTree
+
+def hausdorff_distance(image0, image1):
+    """Calculate the Hausdorff distance between nonzero elements of given images.
+
+    The Hausdorff distance [1]_ is the maximum distance between any point on
+    ``image0`` and its nearest point on ``image1``, and vice-versa.
+
+    Parameters
+    ----------
+    image0, image1 : ndarray
+        Arrays where ``True`` represents a point that is included in a
+        set of points. Both arrays must have the same shape.
+
+    Returns
+    -------
+    distance : float
+        The Hausdorff distance between coordinates of nonzero pixels in
+        ``image0`` and ``image1``, using the Euclidian distance.
+
+    References
+    ----------
+    .. [1] http://en.wikipedia.org/wiki/Hausdorff_distance
+
+    Examples
+    --------
+    >>> points_a = (3, 0)
+    >>> points_b = (6, 0)
+    >>> shape = (7, 1)
+    >>> image_a = np.zeros(shape, dtype=np.bool)
+    >>> image_b = np.zeros(shape, dtype=np.bool)
+    >>> image_a[points_a] = True
+    >>> image_b[points_b] = True
+    >>> hausdorff_distance(image_a, image_b)
+    3.0
+
+    """
+    a_points = np.transpose(np.nonzero(image0))
+    b_points = np.transpose(np.nonzero(image1))
+
+    # Handle empty sets properly:
+    # - if both sets are empty, return zero
+    # - if only one set is empty, return infinity
+    if len(a_points) == 0:
+        return 0 if len(b_points) == 0 else np.inf
+    elif len(b_points) == 0:
+        return np.inf
+
+    return max(max(cKDTree(a_points).query(b_points, k=1)[0]),
+               max(cKDTree(b_points).query(a_points, k=1)[0]))
diff --git a/skimage/metrics/tests/__init__.py b/skimage/metrics/tests/__init__.py
@@ -0,0 +1,9 @@
+from ..._shared.testing import setup_test, teardown_test
+
+
+def setup():
+    setup_test()
+
+
+def teardown():
+    teardown_test()
diff --git a/skimage/metrics/tests/test_set_metrics.py b/skimage/metrics/tests/test_set_metrics.py
@@ -0,0 +1,118 @@
+from __future__ import print_function, division
+
+import numpy as np
+from numpy.testing import assert_almost_equal
+import itertools
+
+from skimage._shared.testing import parametrize
+from skimage.metrics import hausdorff_distance
+
+
+def test_hausdorff_empty():
+    empty = np.zeros((0, 2), dtype=np.bool)
+    non_empty = np.zeros((3, 2), dtype=np.bool)
+    assert hausdorff_distance(empty, non_empty) == 0.
+    assert hausdorff_distance(non_empty, empty) == 0.
+    assert hausdorff_distance(empty, empty) == 0.
+
+
+def test_hausdorff_simple():
+    points_a = (3, 0)
+    points_b = (6, 0)
+    shape = (7, 1)
+    coords_a = np.zeros(shape, dtype=np.bool)
+    coords_b = np.zeros(shape, dtype=np.bool)
+    coords_a[points_a] = True
+    coords_b[points_b] = True
+    distance = np.sqrt(sum((ca - cb) ** 2
+                           for ca, cb in zip(points_a, points_b)))
+    assert_almost_equal(hausdorff_distance(coords_a, coords_b), distance)
+
+
+@parametrize("points_a, points_b",
+             itertools.product([(0, 0), (3, 0), (1, 4), (4, 1)], repeat=2))
+def test_hausdorff_region_single(points_a, points_b):
+    shape = (5, 5)
+    coords_a = np.zeros(shape, dtype=np.bool)
+    coords_b = np.zeros(shape, dtype=np.bool)
+    coords_a[points_a] = True
+    coords_b[points_b] = True
+
+    distance = np.sqrt(sum((ca - cb) ** 2
+                           for ca, cb in zip(points_a, points_b)))
+    assert_almost_equal(hausdorff_distance(coords_a, coords_b), distance)
+
+
+@parametrize("points_a, points_b",
+             itertools.product([(5, 4), (4, 5), (3, 4), (4, 3)],
+                               [(6, 4), (2, 6), (2, 4), (4, 0)]))
+def test_hausdorff_region_different_points(points_a, points_b):
+    shape = (7, 7)
+    coords_a = np.zeros(shape, dtype=np.bool)
+    coords_b = np.zeros(shape, dtype=np.bool)
+    coords_a[points_a] = True
+    coords_b[points_b] = True
+
+    distance = np.sqrt(sum((ca - cb) ** 2
+                           for ca, cb in zip(points_a, points_b)))
+    assert_almost_equal(hausdorff_distance(coords_a, coords_b), distance)
+
+
+def test_gallery():
+    shape = (60, 60)
+
+    # Create a diamond-like shape where the four corners form the 1st set
+    # of points
+    x_diamond = 30
+    y_diamond = 30
+    r = 10
+
+    plt_x = [0, 1, 0, -1]
+    plt_y = [1, 0, -1, 0]
+
+    set_ax = [(x_diamond + r * x) for x in plt_x]
+    set_ay = [(y_diamond + r * y) for y in plt_y]
+
+    # Create a kite-like shape where the four corners form the 2nd set of
+    # points
+    x_kite = 30
+    y_kite = 30
+    x_r = 15
+    y_r = 20
+
+    set_bx = [(x_kite + x_r * x) for x in plt_x]
+    set_by = [(y_kite + y_r * y) for y in plt_y]
+
+    # Set up the data to compute the hausdorff distance
+    coords_a = np.zeros(shape, dtype=np.bool)
+    coords_b = np.zeros(shape, dtype=np.bool)
+
+    for x, y in zip(set_ax, set_ay):
+        coords_a[(x, y)] = True
+
+    for x, y in zip(set_bx, set_by):
+        coords_b[(x, y)] = True
+
+    # Test the hausdorff function on the coordinates
+    # Should return 10, the distance between the furthest tip of the kite and
+    # its closest point on the diamond, which is the furthest someone can make
+    # you travel to encounter your nearest neighboring point on the other set.
+    assert_almost_equal(hausdorff_distance(coords_a, coords_b), 10.)
+
+
+@parametrize("points_a, points_b",
+             itertools.product([(0, 0, 1), (0, 1, 0), (1, 0, 0)],
+                               [(0, 0, 2), (0, 2, 0), (2, 0, 0)]))
+def test_3d_hausdorff_region(points_a, points_b):
+    hausdorff_distances_list = []
+    shape = (3, 3, 3)
+    coords_a = np.zeros(shape, dtype=np.bool)
+    coords_b = np.zeros(shape, dtype=np.bool)
+    coords_a[points_a] = True
+    coords_b[points_b] = True
+
+    distance = np.sqrt(sum((ca - cb) ** 2
+                           for ca, cb in zip(points_a, points_b)))
+    hausdorff_distance_3d = hausdorff_distance(coords_a, coords_b)
+    assert_almost_equal(hausdorff_distance_3d, distance)
+    hausdorff_distances_list.append(hausdorff_distance_3d)
diff --git a/skimage/setup.py b/skimage/setup.py
@@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None):
     config.add_subpackage('graph')
     config.add_subpackage('io')
     config.add_subpackage('measure')
+    config.add_subpackage('metrics')
     config.add_subpackage('morphology')
     config.add_subpackage('transform')
     config.add_subpackage('util')