Distance transforms of sampled functions

We describe linear-time algorithms for solving a class of problems that involve transforming a cost function on a grid using spatial information. These problems can be viewed as a generalization of classical distance transforms of binary images, where the binary image is replaced by an arbitrary function on a grid. Alternatively they can be viewed in terms of the minimum convolution of two functions, which is an important operation in grayscale morphology. A consequence of our techniques is a simple and fast method for computing the Euclidean distance transform of a binary image. Our algorithms are also applicable to Viterbi decoding, belief propagation, and optimal control.


Introduction
Distance transforms are an important tool in computer vision, image processing and pattern recognition.A distance transform of a binary image specifies the distance from each pixel to the nearest "on" pixel.Distance transforms play a central role in the comparison of binary images, particularly for images resulting from local feature detection techniques such as edge or corner detection.For example, both the Chamfer [6] and Hausdorff [16] matching approaches make use of distance transforms in comparing binary images.Distance transforms are also used to compute the medial axis (boundaries of Voronoi cells) of digital shapes [4].
In this paper we consider a generalization of distance transforms to arbitrary functions on a grid rather than binary-valued ones (i.e., real-valued images rather than binary images).There is a simple intuition underlying this generalization.Binary images are often used to specify the locations of certain features in a picture.Rather than using a binary array to specify the presence or absence of a feature at each pixel, it can be useful to have a real-valued array specifying a cost (or strength) for a feature at each pixel.For such more general situation it is again useful to compute a type of distance transform, which in this case should reflect a combination of distances and feature costs.
Let G be a regular grid and f : G → R a function on the grid.We refer to a function on a grid as a sampled function.Sampled functions are usually defined by arrays.We define the distance transform of f to be a another function D f : G → R where D f (p) = min q∈G d(p, q) + f (q) . (1.1) Here d(p, q) is some measure of the distance (not necessarily a metric) between p and q.Intuitively, for each point p we find a point q that is close to p, and for which f (q) is small.Note that if f has a small value at some location, D f will have small value at that location and any nearby point, where nearness is measured by the distance d(p, q).When d(p, q) is the Euclidean distance, we refer to D f as the Euclidean distance transform (EDT) of f .The definition in equation (1.1) is closely related to the traditional distance transform of a set of points P ⊆ G, which associates to each grid location the distance to the nearest point in P, The grid G and the set of points P are often defined by a binary image, in which case the distance transform D P is a real-valued image of the same size.Many algorithms for computing the distance transform of point sets use the equivalent definition where 1(q) is an indicator function for membership in P: This definition is almost the same as the definition for the distance transform of a sampled function in equation (1.1), except that it uses the indicator function 1(q) rather than an arbitrary function f (q).Efficient algorithms for computing the distance transform of a binary image using the L 1 and L ∞ distances were developed by Rosenfeld and Pfaltz [27].Similar methods described in [5] have been widely used to compute approximations to the EDT of a binary image.These algorithms can be easily adapted to compute the distance transform of a sampled function, as we show in Section 3 for the case of the L 1 distance.
Our main result is a new linear-time algorithm for computing the distance transform of a sampled function when distance is measured by the squared Euclidean distance.This in turn provides a new technique for computing the exact EDT of a binary image, by computing the transform of the corresponding indicator function and taking square roots of the result.There are a number of other algorithms for computing the EDT of a binary image in linear-time (e. g., [18,7,17]); however these methods are quite involved and are not widely used in practice.In contrast, our algorithm is relatively simple, easy to implement, and very fast in practice.In the Conclusions (Section 5) we indicate the range of areas of science and technology to which our algorithm has already been applied.
A sampled function is equivalent to a real-valued image.We use the terminology "distance transform of a sampled function" for two reasons.First, we want to stress that the input is generally some kind of cost function that is being transformed so as to incorporate spatial (or distance) information.Second, there are methods for computing transforms of real-valued images based on minimum distances along paths, where the cost of a path is the sum of values along the path [28].We want to avoid confusion with these methods which compute something quite different from what we consider here.

Minimum convolution
The distance transform of a function is closely related to the min-convolution operation.This operation and its continuous counterpart play an important role in grayscale morphology [21].The min-convolution of two functions f and g is defined as, Just like standard convolution, this operation is commutative and associative, When d(p, q) = g(p − q), the distance transform of f under the distance d is the min-convolution of f and g.The algorithms we describe in this paper can be seen as a algorithms for computing minconvolutions with restricted functions g.Our basic algorithms solve one-dimensional problems.The multi-dimensional case is handled by decomposing it into a sequence of one-dimensional problems.
In the case of the squared Euclidean distance we are computing the min-convolution of a onedimensional function f and a parabola.We solve the problem by noting a relationship to lower envelopes of parabolas.Lower envelopes of n parabolas can be computed in O(n log n) time using the method in [8].That method operates by sorting the parabolas in an appropriate order to be inserted in the lower envelope in amortized constant time.In our case we get an O(n) algorithm because the parabolas have the same shape and are already sorted in an appropriate order.Here n is the size of the domain of f .In the case of the L 1 distance we compute the min-convolution of a one-dimensional function f and a cone using the classical approach of Rosenfeld and Pfaltz [27].
When d is the squared Euclidean and L 1 distance, the corresponding g is convex.As described by Eppstein [9] the problem of computing one-dimensional min-convolutions between arbitrary functions f and convex functions g can be solved in O(n) time using the totally monotone matrix search algorithm of [1].When g is concave the problem can be solved in O(nα(n)) time using the matrix search algorithm of [20], where α is the extremely slowly growing inverse Ackermann function.
Our algorithm for one-dimensional min-convolution with a parabola generalizes to computing the min-convolution of an arbitrary one-dimensional function f and a convex or concave function g.If an intersection point between shifted copies of g can be computed in constant time we get an O(n) method for computing the min-convolution (the resulting algorithm is simpler than the linear-time method based on matrix search, and also applies to the case when g is concave).Otherwise binary search can be used to find an intersection in O(log n) time, yielding an O(n log n) min-convolution algorithm.

Optimization problems
Distance transforms of sampled functions arise in the solution of a number of optimization problems.For instance in the widely used Viterbi algorithm for hidden Markov models [24], in max-product belief propagation [23], in optimal control methods [3] and in resource allocation problems [2].
In these problems we typically have a discrete state space S, a cost b(p) for each state p ∈ S, a transition cost a(p, q) for changing from state p to state q, and a dynamic programming equation For our purposes a detailed understanding of this equation is not as important as observing its form.The minimization in the second term of the right hand side is closely related to the distance transform of a function.If S is a grid, then equation (1.2) can be rewritten in terms of a distance transform Thus algorithms for computing distance transforms of functions apply to certain problems of the form in (1.2).We have used these methods to develop improved algorithms for recognition of articulated (non-rigid) objects [11], for inference using large state-space hidden Markov models [13], and for the solution of low-level vision problems such as stereo, image restoration and optical flow using loopy belief propagation [12].For instance, in the case of a hidden Markov model with n states the standard computation of the Viterbi recurrence takes O(n 2 ) time which is not practical for large values of n, while the computation using our distance transform techniques takes O(n) time.

Classical EDT
An important consequence of our results is a new linear-time algorithm for computing the classical Euclidean distance transform of a binary image.The method works by decomposing the problem into a sequence of one-dimensional problems.This decomposition relies on the solution of the more general problem, involving transforms of sampled functions.Thus, by solving a more general problem for the one-dimensional case we obtain a simple algorithm for the multi-dimensional EDT.While other methods for computing the EDT of binary functions existed before, these methods have not been used in practice due to their practical complexity.In contrast, our algorithm has been used in a wide variety of applications.For example, it was applied to study urban change in [15], motion planning in [25], object recognition in [29] and image analysis PDEs in [22].This illustrates the significance of developing algorithms that are not only asymptotically optimal, but are also easy to implement and have low complexity in practice.

Squared Euclidean distance 2.1 One dimension
Let G = {0, . . ., n − 1} be a one-dimensional grid, and f : G → R be a function on the grid.The squared Euclidean (or quadratic) distance transform of f is given by (2.1) THEORY OF COMPUTING, Volume 8 (2012), pp.415-428 Note that for each point q ∈ G there is a constraint that the distance transform of f be bounded from above by a parabola rooted at (q, f (q)).In fact the distance transform is defined by the lower envelope of these parabolas, as shown in Figure 1.The value of the distance transform at p is simply the height of the lower envelope at that point.
Our algorithm for computing this distance transform has two distinct steps.First we compute the lower envelope of the n parabolas just mentioned.We then fill in the values of D f (p) by checking the height of the lower envelope at each grid location p.This is a unique approach because we start with something defined on a grid (the values of f ), move to a combinatorial structure defined over the whole domain (the lower envelope of the parabolas) and move back to values on the grid by sampling the lower envelope.Pseudocode for the procedure is shown in Algorithm 1.
The main part of the algorithm is the lower envelope computation.Note that any two parabolas defining the distance transform intersect at exactly one point.Simple algebra yields the horizontal position of the intersection between the parabolas coming from grid positions q and r as, If q < r then the parabola coming from q is below the one coming from r to the left of the intersection point s, and above it to the right of s.
We compute the lower envelope by sequentially computing the lower envelope of the first q parabolas, where the parabolas are ordered according to the horizontal locations of their vertices.The algorithm works by computing the combinatorial structure of this lower envelope.We keep track of the structure by using two arrays.The horizontal grid location of the i-th parabola in the lower envelope is stored in v[i].The range in which the i-th parabola of the lower envelope is below the others is given by z[i] and z[i + 1].The variable k keeps track of the number of parabolas in the lower envelope.
When considering the parabola from q, we find its intersection with the parabola from v[k] (the rightmost parabola in the lower envelope computed so far).There are two possible cases, as illustrated in s ← (( Algorithm 1: One-dimensional distance transform under the squared Euclidean distance.
from q is below all others starting at the intersection point.If the intersection is before z[k] then the parabola from v[k] should not be part of the new lower envelope, so we decrease k to delete that parabola and repeat the procedure.
Theorem 2.1.Algorithm 1 correctly computes the distance transform of a one-dimensional sampled function under the squared Euclidean distance in O(n) time.
Proof.We start by showing that the algorithm correctly computes the lower envelope of the first q parabolas by induction.The base case is trivial.The lower envelope of the first parabola is represented by a single interval from −∞ to +∞ dominated by the parabola from grid position 0. Let s be the horizontal position of the intersection between the q-th parabola and the one from v[k] as computed in line 6.The parabola from q is above the one from v[k] to the left of s, and below it to the right of s.By the induction hypothesis the parabola from v[k] is above at least one other parabola in the lower envelope to the left of z[k] and below all of them to the right of z[k].
Suppose s > z[k] (as in Figure 2a).The parabola from q is below the one from v[k] to the right of s, which in turn is below all others everywhere after z[k].We see that the parabola from q dominates the lower envelope after s.The parabola from v[k] continues to be the lowest between z[k] and s.The algorithm updates the lower envelope to reflect these changes by creating a new interval from s to ∞ dominated by the q-th parabola.
Suppose s ≤ z[k] (as in Figure 2b).Now the parabola from q is below the parabola from v[k] in the whole interval previously dominated by the parabola from v[k].This means that the parabola from v

Figure 2:
The two possible cases considered by the algorithm when adding the parabola from q to the lower envelope constructed so far.
is not part of the lower envelope of the first q parabolas.The algorithm modifies the lower envelope to remove the parabola from v[k] and proceeds to add the parabola from q to the remaining structure.This process eventually terminates since z Once the lower envelope is computed it remains to fill in the distance transform values by sampling the height of the lower envelope at each grid location.This is done from left to right on lines 14 through 18.To understand the running time of the algorithm note that we consider adding each parabola to the lower envelope exactly once.The addition of one parabola may involve the deletion of many others, but each parabola is deleted at most once.So the overall computation of the lower envelope in lines 1 through 13 takes O(n) time.The computation of the distance transform values from the lower envelope in lines 14 through 18 considers each grid position and each parabola in the lower envelope at most once, so the second part of the algorithm also takes O(n) time.

Arbitrary dimensions
Let G = {0, . . ., n − 1} × {0, . . ., m − 1} be a two-dimensional grid, and f : G → R a function on the grid.The two-dimensional distance transform of f under the squared Euclidean distance is given by, The first term does not depend on y so we can rewrite this equation as, Thus a two-dimensional distance transform can be computed by first computing one-dimensional distance transforms along each column of the grid, and then computing one-dimensional distance transforms along each row of the result.
This argument extends to arbitrary dimensions, resulting in the composition of transforms along each dimension of the underlying grid.Note that changing the order of these transforms yields the same result, as can be seen readily for the two-dimensional case above.The multi-dimensional algorithm runs in O(dN) time, where d is the dimension of the grid and N is the overall number of grid locations (d = 2 and N = nm for the grid defined above).Therefore we obtain the following result.
Theorem 2.2.The squared Euclidean distance transform of a sampled function defined on a d dimensional grid with N sample points can be computed in O(dN) time.
Figure 3 illustrates the computation of the classical EDT of a binary image using our method, where we start with the indicator function for the points on the grid.Note that in equation (2.3) the function that must be transformed along the second dimension is no longer an indicator function.Thus the notion of a distance transform for arbitrary functions is important here.The separation of the multi-dimensional transform into multiple one-dimensional ones makes our method substantially simpler than previous techniques for computing the classical EDT.

L 1 distance
For a set of points on a one-dimensional grid, the distance transform under the L 1 distance can be computed using a simple two-pass algorithm (e. g., [27]).Essentially the same algorithm can be used to compute the distance transform of a one-dimensional sampled function under the L 1 distance.Pseudocode for the method is shown in Algorithm 2.
The values of the distance transform are initialized to the values of f itself.In the forward pass, each successive element of D f (q) is set to the minimum of its own value and one plus the value of the previous element.This is done "in place" so that updates affect one another.In the backward pass each item is analogously set to the minimum of its own value and one plus the value of the next element.For example given the input (4, 2, 8, 6, 1), after the first pass D f will be (4, 2, 3, 4, 1), and after the second pass it will be (3, 2, 3, 2, 1).
[recursive substitution] The last equation is exactly the computation performed by the algorithm.The convolution of the resulting function with b in the backward pass is analogous.Both the forward and backward passes take a constant number of operations per grid position, yielding a O(n) algorithm overall.
As with the squared Euclidean distance case considered in the previous section, a multi-dimensional L 1 distance transform can be computed by successive transforms along each dimension of the grid.Theorem 3.2.The L 1 distance transform of a sampled function defined on a d dimensional grid with N sample points can be computed in O(dN) time.

Other distances
We can compute distance transforms of functions under other distances not discussed so far.An important case is the distance transform under the box distance defined by d(p, q) = 0 when |p − q| < T and ∞ otherwise (here T is a threshold).This transform can be computed using a linear time algorithm for the min-filter described in [14] but our approach is simpler and more general.
Algorithm 1 can be easily modified to use any distance of the form d(p, q) = g(p − q) as long as g is convex.We only need to be able to compute intersection points between shifted copies of g defined by g q (p) = f (q) + g(p − q) for q ∈ {0, . . ., n − 1}.This can often be done in constant time if g is given in analytic form.In the worst case an intersection point can be found in O(log n) time using binary search.
The key is to note that D f (p) = min q g q (p) and for any convex function g, if q < r we have that g q (p) is below g r (p) before their intersection point, and above afterwards.
Algorithm 1 can also be modified to handle the case when g is concave.The only difference is that we have to consider the shifted copies of g in reverse order.This is because when g is concave and q < r we have that g q (p) is above g r (p) up to their intersection point and below afterwards.
These observations lead to the following result in terms of min-convolutions.
Theorem 4.1.Let f be an arbitrary one-dimensional function on a grid of size n, and g be a concave or convex function.The min-convolution f ⊗ g can be computed in O(nt(n)) time where t(n) is the time required to computed an intersection between shifted copies of g.
Another way to obtain fast distance transform (and min-convolution) algorithms is to use the relationships described below.For example, the distance d(p, q) = min c(p − q) 2 , a|p − q| + b is commonly used in robust estimation and is very important in practice.Intuitively this distance is robust because it grows slowly after a certain point.We can compute the distance transform of a function under the robust distance by computing both a linear and a quadratic distance transform and combining the results.Many state transition costs which appear in dynamic programming equations similar to (1.2) can be written as a combination of a small number of linear and quadratic and box distances using these relations.In such cases our algorithms can be used to compute the dynamic programming equation in time linear in the number of possible states.

Conclusions
Classical Euclidean distance transforms of binary images have many important applications.Previous exact algorithms have not been used in practice because they were too complex despite having linear-time THEORY OF COMPUTING, Volume 8 (2012), pp.415-428 worst case running time.In contrast we have described simple methods for solving more general problems.These problems are equivalent to min-convolutions with special functions.The results described here first appeared in a technical report [10].
Our approach has lead to a practical algorithm for the classical EDT that has been widely used by other researchers.Our EDT algorithm has been used in robotics for motion planning [25], in cartography for understanding urban development [15], in image processing for image segmentation [22], in medical image analysis for prenatal genetic testing [26], and in information interfaces for sonification of images for the visually impared [30].
The generalization of the distance transform to real-valued functions has proven to have important practical applications as well.It has been used in computer vision for object recognition [11] and depth estimation [12] as well as in machine learning for maximum a posteriori inference [13,19].

Figure 1 :
Figure 1: The distance transform as the lower envelope of n parabolas.

Figure 3 :
Figure 3: Illustration of the distance transform of a binary image using the squared Euclidean distance.The input image is shown in (a).The transform of the input along each column is shown in (b).The final distance transform, obtained by transforming the intermediate result along each row, is shown in (c).Dark pixels correspond to low values while bright pixels correspond to high values.

Theorem 3 . 1 .
Algorithm 2 correctly computes the distance transform of a one-dimensional sampled function under the L 1 distance in O(n) time.Proof.Let a(p) = |p| if p ≥ 0 , ∞ otherwise, and b(p) = |p| if p ≤ 0 , ∞ otherwise.It's not hard to check that (a ⊗ b)(p) = |p|.We claim that the forward pass of the algorithm computes the minimum convolution of f with a while the backward pass computes the minimum convolution of the result with b.Since minimum convolution is an associative operation the algorithm computes f ⊗ (a ⊗ b), which as discussed in Section 1.1 is the distance transform of f under the L 1 distance.Now we show that the forward pass of the algorithm does indeed compute the minimum convolution of f and a.