Fastest Algorithm for Distance Transform

Question

I am looking for the fastest available algorithm for distance transform.

According to Image Processing Learning Resources - HIPR 2 (HYPERMEDIA IMAGE PROCESSING REFERENCE) - Morphology - Distance Transform:

The distance transform can be calculated much more efficiently using clever algorithms in only two passes (e.g. Rosenfeld and Pfaltz 1968).

Searching around, I found: "Rosenfeld, A and Pfaltz, J L. 1968. Distance Functions on Digital Pictures. Pattern Recognition, 1, 33-61."

But I believe we should have a better and faster algorithm than the one in 1968 already? In fact, I could not find the source from 1968, so any help is highly appreciated.

Answer 1

Pedro F. Felzenszwalb and Daniel P. Huttenlocher have published their implementation for the distance transform [archive]. You cannot use it for volumetric images, but maybe you can extend it to support 3d data. I have only used it as a black box.

Answer 2

This paper discusses all of the modern exact distance transforms:

"2D Euclidean distance transforms: a comparative survey", ACM Computing Surveys, Vol 40, Issue 1, Feb 2008 http://www.lems.brown.edu/~rfabbri/stuff/fabbri-EDT-survey-ACMCSurvFeb2008.pdf

The paper cites the technique from Meijster, et. al. as the fastest general purpose, exact transform. This technique is detailed here:

"A General Algorithm for Computing Distance Transforms in Linear Time", A. Meijster, J. B. T. M. Roerdink and W. H. Hesselink. http://fab.cba.mit.edu/classes/S62.12/docs/Meijster_distance.pdf

The Meijster algorithm is used in my open source effects library: https://github.com/vinniefalco/LayerEffects

I hope this helps someone.

Answer 3

Here is a C# code for 1D squared euclidean distance transform according to the Felzenszwald & Huttenlocher's paper:

private static void DistanceTransform(double[] dataInput, ref double[] dataOutput)
{
    int n = dataInput.Length;

    int k = 0;
    int[] v = new int[n];
    double[] z = new double[n + 1];

    v[0] = 0;
    z[0] = Double.NegativeInfinity;
    z[1] = Double.PositiveInfinity;

    double s;

    for (int q = 1; q < n; q++)
    {
        while (true)
        {
            s = (((dataInput[q] + q * q) - (dataInput[v[k]] + v[k] * v[k])) / (2.0 * q - 2.0 * v[k]));

            if (s <= z[k])
            {
                k--;
            }
            else
            {
                break;
            }
        }

        k++;

        v[k] = q;
        z[k] = s;
        z[k + 1] = Double.PositiveInfinity;
    }

    k = 0;

    for (int q = 0; q < n; q++)
    {
        while (z[k + 1] < q)
        {
            k++;
        }

        dataOutput[q] = ((q - v[k]) * (q - v[k]) + dataInput[v[k]]);
    }
}

This can be readily used for binary and grayscale images by applying it first on image columns and then rows (or vice versa, of course).

The transform is indeed very fast.

Here are the source and output images:

The black pixels have value 0 and the white have some large value (have to be larger than largest possible squared distance in the images but not infinity) so that the transform returns distance from the black pixels and the white ones are ommited.

To get true euclidean distance transform, simply take a square root of each pixel from the output image.