7

I have multiple grids (numpy arrays [Nk,Ny,Nx]) and would like to use Hausdorff distance as a metric of similarity of these grids. There are several modules in scipy (scipy.spatial.distance.cdist,scipy.spatial.distance.pdist) which allow to calculate Euclidean distance between 2D arrays. Now to compare grids I have to choose some cross-section (e.g. grid1[0,:] & grid2[0,:]) and compare it between each other. Is it possible to calculate Hausdorff distance between 3D grids directly?

Vitali Molchan
  • 203
  • 4
  • 9
  • This question (http://stackoverflow.com/questions/13692801/distance-matrix-of-curves-in-python) may be relevant. The conclusion seem to be there is no scipy/numpy algorithm and that it is better to write the main algorithm in c if speed is critical (that was for 2D). – Ed Smith Jun 08 '15 at 11:37
  • Farhawa, thanks a lot! – Vitali Molchan Jun 08 '15 at 14:06

3 Answers3

8

I am newby here, but faced with the same challenge and tried to attack it directly on a 3D level.

So here is the function I did:

def Hausdorff_dist(vol_a,vol_b):
dist_lst = []
for idx in range(len(vol_a)):
    dist_min = 1000.0        
    for idx2 in range(len(vol_b)):
        dist= np.linalg.norm(vol_a[idx]-vol_b[idx2])
        if dist_min > dist:
            dist_min = dist
    dist_lst.append(dist_min)
return np.max(dist_lst)

The input needs to be numpy.array, but the rest is working directly.

I have 8000 vs. 5000 3D points and this runs for several minutes, but at the end it gets to the distance you are looking for.

This is however checking the distance between two points, not neccesarily the distance of two curves. (neither mesh).

Edit (on 26/11/2015):

Recenty finished the fine-tuned version of this code. Now it is splitted into two part.

First is taking care of grabbing a box around a given point and taking all the radius. I consider this as a smart way to reduce the number of points required to check.

def bbox(array, point, radius):
    a = array[np.where(np.logical_and(array[:, 0] >= point[0] - radius, array[:, 0] <= point[0] + radius))]
    b = a[np.where(np.logical_and(a[:, 1] >= point[1] - radius, a[:, 1] <= point[1] + radius))]
    c = b[np.where(np.logical_and(b[:, 2] >= point[2] - radius, b[:, 2] <= point[2] + radius))]
    return c

And the other code for the distance calculation:

def hausdorff(surface_a, surface_b):

    # Taking two arrays as input file, the function is searching for the Hausdorff distane of "surface_a" to "surface_b"
    dists = []

    l = len(surface_a)

    for i in xrange(l):

        # walking through all the points of surface_a
        dist_min = 1000.0
        radius = 0
        b_mod = np.empty(shape=(0, 0, 0))

        # increasing the cube size around the point until the cube contains at least 1 point
        while b_mod.shape[0] == 0:
            b_mod = bbox(surface_b, surface_a[i], radius)
            radius += 1

        # to avoid getting false result (point is close to the edge, but along an axis another one is closer),
        # increasing the size of the cube
        b_mod = bbox(surface_b, surface_a[i], radius * math.sqrt(3))

        for j in range(len(b_mod)):
            # walking through the small number of points to find the minimum distance
            dist = np.linalg.norm(surface_a[i] - b_mod[j])
            if dist_min > dist:
                dist_min = dist

        dists.append(dist_min)

    return np.max(dists)
Akos Gulyban
  • 121
  • 1
  • 8
1

In case anyone is still looking for the answer to this question years later... since 2016 scipy now includes a function to calculate the Hausdorff distance in 3D:

scipy.spatial.distance.directed_hausdorff

craq
  • 1,441
  • 2
  • 20
  • 39
  • this function is 2D only... – Guillaume Mougeot Jan 06 '23 at 10:54
  • interesting... the documentation says that it is 2D only, but it works in 3D for my tests. Looking at [the code](https://github.com/scipy/scipy/blob/dde50595862a4f9cede24b5d1c86935c30f1f88a/scipy/spatial/_hausdorff.pyx) there doesn't seem to be any limitation on array dimensions at all. Do you agree? If so, I might raise a ticket with scipy for them to update their documentation. – craq Jan 08 '23 at 23:27
  • 1
    you might actually seem to be right! Yes, maybe raise a ticket with scipy – Guillaume Mougeot Jan 09 '23 at 17:47
  • No need to have a ticket. The [function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.directed_hausdorff.html) returns the “distance between two 2-D arrays of coordinates”, where 2-D means a 2D array. You may have 2 lists of N-dimensional points. – Dom May 11 '23 at 15:49
1

There is the package point_cloud_utils which provides a few 3D metrics, such as Hausdorff distance. Install it with pip install point_cloud_utils and then use it like follows:

import point_cloud_utils as pcu
import numpy as np

# Generate two random point sets
a = np.random.rand(1000, 3)
b = np.random.rand(500, 3)

# Compute one-sided squared Hausdorff distances
hausdorff_a_to_b = pcu.one_sided_hausdorff_distance(a, b)
hausdorff_b_to_a = pcu.one_sided_hausdorff_distance(b, a)

# Take a max of the one sided squared  distances to get the two sided Hausdorff distance
hausdorff_dist = pcu.hausdorff_distance(a, b)
Philipp
  • 652
  • 2
  • 10
  • 28