University of Groningen
A general algorithm for computing distance transforms in linear time
Meijster, A.; Roerdink, J.B.T.M.; Hesselink, W.H.
Published in:
MATHEMATICAL MORPHOLOGY AND ITS APPLICATIONS TO IMAGE AND SIGNAL PROCESSING
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from
it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date:
2000
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Meijster, A., Roerdink, J. B. T. M., & Hesselink, W. H. (2000). A general algorithm for computing distance
transforms in linear time. In J. Goutsias, L. Vincent, & DS. Bloomberg (Eds.),
MATHEMATICAL
MORPHOLOGY AND ITS APPLICATIONS TO IMAGE AND SIGNAL PROCESSING
(pp. 331-340).
(COMPUTATIONAL IMAGING AND VISION; Vol. 18). University of Groningen, Johann Bernoulli Institute
for Mathematics and Computer Science.
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the
author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.
More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-
amendment.
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the
number of authors shown on this cover page is limited to 10 maximum.
Download date: 10-08-2022
A GENERAL ALGORITHM FOR COMPUTING DISTANCE
TRANSFORMS IN LINEAR TIME
A. MEIJSTER
*
‚ J.B.T.M. ROERDINK
†
and W.H. HESSELINK
†
University of Groningen
P.O. Box 800, 9700 AV Groningen, The Netherlands
email: a. meijster@rc.rug. nl,
{
roe, wim
}
@cs.rug.nl
Abstract. A new general algorithm for computing distance transforms of digital images
is presented. The algorithm consists of two phases. Both phases consist of two scans, a
forward and a backward scan. The first phase scans the image column-wise, while the second
phase scans the image row-wise. Since the computation per row (column) is independent of
the computation of other rows (columns), the algorithm can be easily parallelized on shared
memory computers. The algorithm can be used for the computation of the exact Euclidean,
Manhattan (L
1
norm), and chessboard distance (L
∞
norm) transforms.
Key words: Distance Transforms, Row-Column Factorization, Parallelization.
1. Introduction
Distance transforms play an important role in many morphological image proce-
ssing applications. They have been extensively studied and used in computa-
tional geometry, image processing, computer graphics and pattern recognition,
e.g., [1, 2, 3, 7]. The two-dimensional distance transform can be described as
follows. Let B be a set of grid points taken from a rectangular grid of size
m × n . The problem is to assign to every grid point (x , y) the distance to the
nearest point in B. If we use the Euclidean metric for computing distances,
and represent B by a boolean array b[·, ·], we thus want to compute the two
dimensional array dt[x, y ] = , where
Here we use the notation for the minimal value of f( k )
when k ranges over all values that satisfy P(k).
Since the exact Euclidean distance transform is often regarded as too com-
putationally intensive, several algorithms have been proposed that use some
mask which is swept over the image in two scans, to compute approximations
like the Manhattan (city-block) distance, the chessboard distance, or chamfer
distances (see [1, 2, 3, 7]). The time complexity is linear in the number of
pixels of the image (i.e. O(m × n )), but it does not yield the exact Euclidean
distance, which is required for some applications. Another drawback of these
*
A. Meijster works at the Computing Centre of the University of Groningen.
†
J.B.T.M. Roerdink and W.H. Hesselink work at the Institute for Mathematics and Com-
puting Science.
332
A. MEIJSTER ET AL.
algorithms is that they are hard to parallelize for parallel computers since pre-
viously computed results are propagated during the computation, making the
process highly sequential. A recursive algorithm of order mn log m for the ex-
act EDT is given in [5]. In [6] a recursive algorithm of order mn for the exact
EDT is given by reducing the problem to a matrix search algorithm.
In this paper, which is based upon [4], we present a new algorithm that
also computes distance transforms in linear time, is simpler and more efficient
than [6], and is easy to parallelize. It can compute the Euclidean (EDT), the
Manhattan (MDT), and the chessboard distance (CDT) transform, defined by
If we define the minimum of the empty set to be
∞, and use the rule z + ∞ = ∞
for all z, we find with some calculation
where
The algorithm can be summarized as follows. In a first phase each column
C
x
(defined by points (x, y) with x fixed) is separately scanned. For each point
(
x
, y
) on
C
x
, the distance G(x, y) of (x, y ) to the nearest points of C
x
∩ B is
determined. In a second phase each row R
y
(defined by points (x , y ) with y
fixed) is separately scanned, and for each point (x, y ) on R
y
the minimum of
for EDT, for MDT, and
max
for CDT is determined, where (x', y ) ranges over row R
y
.
2. The First Phase
The object of the first phase is to determine the function G. We first observe
that we can split G into two functions GT (top) and GB (bottom), such that
min
where
We start with the computation of GT by introducing an array g to store its
values. It is easy to see that GT(i, y ) = 0 if b [i , y ] holds, and that, otherwise,
GT(
i, y ) = GT(
i
, y – 1) + 1 (or ∞ if y = 0). We can therefore compute
g[ x, y
] := GT(
x, y) using only g[x, y – 1]
in a simple column scan from top to
bottom. Similarly, we find GB(i, y ) = GB(i‚ y + 1) + 1. The second scan runs
from bottom to top, and computes G(x, y ) directly, using GT from the previous
scan, and GB from the current one. After some simplification, this results in
the code fragment given in Fig. 1. Clearly, the time complexity is linear in the
A GENERAL ALGORITHM FOR COMPUTING DISTANCE TRANSFORMS
333
Fig. 1.
Program fragments for the first phase.
number of pixels (i.e. O
(
m × n )). In actual implementations it is convenient
to replace
∞ by m + n, since all distances in the images are less than m + n if
the set B is non-empty.
3. The Second Phase
In the second phase we want to compute EDT, MDT, or CDT row by row,
i.e. for all x with fixed y. Therefore, in this section we regard y as a constant
and omit it as a parameter in auxiliary functions, and introduce g( i ) = G ( i , y ).
Instead of developing an algorithm for each metric separately, we aim at a more
general algorithm for
(1)
The choice of the function f depends on the metric we wish to use, i.e.
It is helpful to introduce a geometrical interpretation of the minimization
problem of Eq. (1). For any i with 0
≤ i < m , denote by F
i
the function
on the real interval [0, m – 1]. We call i the index of F
i
. In the case
of EDT, the graph of F
i
is a parabola with vertex at (i‚ g (i)). In the case of
334
A. MEIJSTER ET AL.
(a) EDT
(b) MDT (c) CDT
Fig. 2. DT as the lower envelope (solid line) of curves F
i
, 0
≤
i < m (dotted lines). The
dashed vertical lines indicate the transitions between regions.
MDT the parabolas are replaced by V-shaped approximations, while in the case
of CDT we deal with ‘topped off’ V-shaped approximations (see Fig. 2). We
can interpret DT geometrically as the lower envelope of the collection {F
i
⏐0 ≤
i < m } evaluated at integer coordinates, cf. Fig. 2. The lower envelopes
consist of a number of consecutive curve segments, whose index we denote by
s
[0]‚ s
[1]‚...,
s
[
q
] counting from left to right. The projections of the segments
on the x-axis are called regions , and form a partition of the interval [0, m) by
consecutive segments. The computation of DT now consists of two scans. In a
forward (left-to-right) scan the set of regions is determined using an incremental
algorithm. In a backward (right-to-left) scan the values DT(x , y ) are trivially
computed for all x.
We start by replacing the upper bound m in (1) by a variable u and define
The geometric interpretation is that we restrict the set B to the half plane to
the left of u. Clearly, DT(x, y ) = FL(x, m
).
For given upper bound u > 0, we define an index h to be a minimizer at
x if, in the expression for FL(x, u ), the minimal value of f(x, i ) occurs at h.
In general, x may have more than one minimizer. defined as the least index h
with 0
≤ h < u such that f (x, h) ≤ f ( x , i) for all i in the same range, i.e.
(2)
We clearly have FL(x, u ) = f ( x, H( x , u )), hence DT(x , y) = f ( x
,
H
( x
,
m
)).
Therefore, the problem reduces to the computation of H(x , m ).
We consider the sets S
(
u) of the least minimizers that occur during the scan
from left to right, and the sets T (h, u ) of points with the same least minimizer
h
. We thus define
(3)
