Hausdorff distance between 3D grids

Question

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?

Answer 1

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)