Euclidean distance map of mesh (distance to surface)

mesh
1

Hello everyone,

I am looking for a method to calculate the Euclidean Distance Map (EDM) of a mesh, also called Euclidean distance transform.
EDM(x) of a point x within the domain is simply the shortest distance from this point to the domain’s boundary.
Replace ‘point’ with ‘vertece’ (or ‘DOF’ - but my function order is 1, so vertices would suffice) and ‘domain’ with ‘mesh’ and that’s what i would like to calculate.

The image joined is an illustration of a EDM.

Ideally, I am looking for a method that works in parallel (mpi) but it you have a hint about how to do it in serial, that would help me to start with.
I am working with meshes that contain dozens of millions of verteces, so a scalable method is more than welcome:)

Best,

Distance_map

2

Hi, one option (which works in parallel out of the box) is to get the distance by solving the Eikonal equation, see Johan Hake’s answer here. However, if the boundary is complex there are some difficulties in getting the Newton solver to converge on the nonlinear problem that is the Eikonal equation. You can try to make it more robust following the discussion in the linked post. An alternative presented below is to decompose the boundary, solve simpler problems and then get the final distance as the minimum of the distance fields from the boundary pieces.

from dolfin import *


def distance_from_bdry_piece(facet_f, tag):
    '''Solve Eikonal equation to get distance from tagged facets'''
    mesh = facet_f.mesh()
    V = FunctionSpace(mesh, 'CG', 1)
    bc = DirichletBC(V, Constant(0.0), facet_f, tag)

    u, v = TrialFunction(V), TestFunction(V)
    phi = Function(V)
    f = Constant(1.0)

    # Smooth initial guess
    a = inner(grad(u), grad(v))*dx
    L = inner(f, v)*dx
    solve(a == L, phi, bc)

    # Eikonal equation with stabilization
    eps = Constant(mesh.hmax()/100)
    F = sqrt(inner(grad(phi), grad(phi)))*v*dx - inner(f, v)*dx + eps*inner(grad(phi), grad(phi))*v*dx
    solve(F == 0, phi, bc)

    return phi


def function_min(foos):
    '''Function that is min of foos'''
    f = foos.pop()
    V = f.function_space()

    vec = lambda u: as_backend_type(u.vector()).vec()
    # The work vector
    min_f = f.copy()
    f_vec = vec(min_f)
    # Reduce
    while foos:
        g = foos.pop()
        f_vec.pointwiseMin(f_vec, vec(g))
    # Sync
    as_backend_type(min_f.vector()).update_ghost_values()
        
    return min_f


mesh = UnitSquareMesh(32, 32)
facet_f = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
CompiledSubDomain('near(x[0], 0)').mark(facet_f, 1)
CompiledSubDomain('near(x[0], 1)').mark(facet_f, 2)
CompiledSubDomain('near(x[1], 0)').mark(facet_f, 3)
CompiledSubDomain('near(x[1], 1)').mark(facet_f, 4)

phis = [distance_from_bdry_piece(facet_f, tag) for tag in (1, 2, 3, 4)]
for i, phi in enumerate(phis):
    File('d%d.pvd' % i) << phi
# Final distance is the min of these
phi = function_min(phis)
File('distance.pvd') << phi
1 Like