## The distance transform, erosion and separability

David Coeurjolly from the Université de Lyon just gave a presentation here at my department. He discussed, among other things, an algorithm that computes the Euclidean distance transform and is separable. The distance transform is an operation that takes a binary image as input, and writes in each object pixel the distance to the nearest background pixel. All sorts of approximations exist, using various distance measures that approximate the Euclidean distance. Using truly Euclidean distances is rather expensive. However, by making an algorithm that is separable, the computational cost is greatly reduced.

David’s algorithm computes the square distance first along each of the rows of the image, then modifies these distances by doing some operations along the columns. In higher-dimensional images you can just repeat this last step along the other dimensions. The operation to modify these distance values sounded very much like parabolic erosions to me, so I just gave this a try.

First, we create some test image and compute its distance transform. To execute these commands in MATLAB you need DIPimage.

a = threshold(gaussf(readim,5)) b = dt(a,1,'true') dipmapping('lin')

By specifying `'true'`

as an option to the `dt`

function we make sure there’s no errors introduced by speeding up the computation.

Now, to compute the distance transform using the parabolic erosion, we will generate an image that is `Inf`

everywhere inside the object. Then we compute the erosion, which will reduce these initial values according to a quadratic distance to the region with zero values. Note that the erosion with a parabolic structuring element is separable.

c = newim(a); c(a) = Inf; c = erosion(c,1,'parabolic'); c = sqrt(c) dipmapping('lin')

Let’s compare these two results.

any(b~=c)

ans = 0

So this is an exact result. Is it really faster? To time these operations I use Steve Eddins’ `timeit`

function. You can get it from MATLAB Central File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/18798. Also, this time I’m using an input image with a really large values instead of `Inf`

. This is just so I can write the operation as a single MATLAB command: `a*1e9`

produces a valid input image, but `a*Inf`

produces `NaN`

values where `a`

is 0. As long as the value used is larger than the square of the maximal distance in the image, this will not affect the output.

f = @() dt(a,1,'true'); g = @() sqrt(erosion(a*1e9,1,'parabolic')); timeit(f) timeit(g)

ans = 0.0253 ans = 0.0112

Yes, it is quite a bit faster. But of course, there are approximations that are much faster, and not bad either.

f = @() dt(a,1,'ties'); g = @() dt(a,1,'fast'); timeit(f) timeit(g)

ans = 0.0062 ans = 0.0040

As I wrote this, I decided to check with Rein van den Boomgaard’s PhD thesis (written in 1992). Rein describes parabolic morphology as the morphological equivalent of the Gaussian convolution kernel, because it is the only structuring element that is rotationally invariant and separable. He says that Sternberg (1980) showed that the distance transform can be calculated exactly by erosion with a cone structuring element. He then concludes that you can compute the squared of the distance transform using the parabolic structuring element. This means that this separable algorithm to compute the exact Euclidean distance transform has been around since at least 1992. And I should have known about it!

I think the real advantage of this technique is that a separable algorithm can be split into many threads, and thus obtain a considerable speedup on multi-processor and multi-core computers, as discussed by David today.

Wow was i surprised that my thesis is still being read…

Somehow the observation in my thesis that Euclidean DT’s can be calculated with a (separable) parabolic erosion had gone unnoticed for some time. And given the work of Sternberg it was not that important too… Maybe your argument about separability and multi processor is the reason for this (parallel processors were not that important then). The same idea (with very nice algorithms) can be found in the papers of Roerdink en Meijster from Groningen University, they were looking for parallel algorithms!

The problem of course with these separable transforms is that you are not able to calculate the constrained (geodesic) DT with this (as far as i am aware of).

Hi Cris,

Nice analysis and thanks for the reference on Rein’s work…. I didn’t know his work..

However, I’d like to correct some points on the literature: I definitely agree with Rein, computing the DT with a linear time and separable process is a solved problem. The first separable linear in time algorithm have been proposed by Breu, H. and Gil, J. and Kirkpatrick, D. and Werman, M. (Linear Time Euclidean Distance Transform Algorithms, IEEE TPAMI, 1995). Before Roerdnik’s technique, T. Hirata (A unified linear-time algorithm for computing distance maps, Information Processing Letters, 1996) has also proposed a separable process dealing with the square of the DT (Breu’s algorithm used a discrete Voronoi diagram computation). Independently, Roerdink en Meijster (2000) rediscovered Hirata’s result.

To point of my presentation was to demonstrate that, within the same framework (and with the same computational cost), we can also extract a discrete medial axis and perform distance reconstruction in linear time (plus some advantages on parallel architecture and toric domains).

David

Hi David,

Thanks for pointing out those references, I’ll enjoy reading them.

I hope I didn’t misrepresent your presentation. I got excited about the one detail that tied in so neatly with some stuff I had been working on earlier, and wanted to explore that. I never intended to give a full overview of your talk. And, I’m happy to say, separable DTs is not the only thing I learned from it!

Cris.