Hesselink, W.H.

Meijster, A.

Roerdink, J.B.T.M.

English

A new general algorithm fur 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-infinity norm) transforms.</p

ARTS repository - University of Groningen

   University of GroningenA general algorithm for computing distance transforms in linear timeMeijster, A.; Roerdink, J.B.T.M.; Hesselink, W.H.Published in:MATHEMATICAL MORPHOLOGY AND ITS APPLICATIONS TO IMAGE AND SIGNAL PROCESSINGIMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite fromit. Please check the document version below.Document VersionPublisher's PDF, also known as Version of recordPublication date:2000Link to publication in University of Groningen/UMCG research databaseCitation for published version (APA):Meijster, A., Roerdink, J. B. T. M., & Hesselink, W. H. (2000). A general algorithm for computing distancetransforms in linear time. In J. Goutsias, L. Vincent, & DS. Bloomberg (Eds.), MATHEMATICALMORPHOLOGY AND ITS APPLICATIONS TO IMAGE AND SIGNAL PROCESSING (pp. 331-340).(COMPUTATIONAL IMAGING AND VISION; Vol. 18). University of Groningen, Johann Bernoulli Institutefor Mathematics and Computer Science.CopyrightOther 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 theauthor(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 policyIf you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediatelyand investigate your claim.Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons thenumber of authors shown on this cover page is limited to 10 maximum.Download date: 06-06-2022A GENERAL ALGORITHM FOR COMPUTING DISTANCETRANSFORMS IN LINEAR TIMEA. MEIJSTER*‚  J.B.T.M.  ROERDINK† and W.H. HESSELINK†University of GroningenP.O. Box 800, 9700 AV Groningen, The Netherlandsemail: a. meijster@rc.rug. nl, {roe, wim} @cs.rug.nlAbstract. A new general algorithm for computing distance transforms of digital imagesis presented. The algorithm consists of two phases. Both phases consist of two scans, aforward and a backward scan. The first phase scans the image column-wise, while the secondphase scans the image row-wise. Since the computation per row (column) is independent ofthe computation of other rows (columns), the algorithm can be easily parallelized on sharedmemory 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. IntroductionDistance 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 asfollows. Let B be a set of grid points taken from a rectangular grid of sizem × n . The problem is to assign to every grid point (x , y) the distance to thenearest 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 twodimensional array dt[x , y ] = , whereHere 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 somemask which is swept over the image in two scans, to compute approximationslike the Manhattan (city-block) distance, the chessboard distance, or chamferdistances (see [1, 2, 3, 7]). The time complexity is linear in the number ofpixels of the image (i.e. O (m × n )), but it does not yield the exact Euclideandistance, 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 theprocess  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 exactEDT is given by reducing the problem to a matrix search algorithm.In this paper, which is based upon [4], we present a new algorithm thatalso computes  distance  transforms in linear time, is simpler and more efficientthan [6], and is easy to parallelize. It can compute the Euclidean (EDT), theManhattan (MDT), and the chessboard  distance (CDT) transform,  defined byIf we define the minimum of the  empty set to be ∞ , and use the rule z + ∞ = ∞for all z , we find with some calculationwhereThe algorithm  can be summarized as follows. In a first phase  each columnC x (defined  by points  (x , y) with x fixed) is separately scanned. For each point(x, y) on Cx , the distance G(x , y ) of (x , y ) to the nearest points of C x ∩ B i sdetermined. In a second phase each row R y (defined by points (x , y ) with yfixed) is separately scanned, and for each point (x , y ) on R y the minimum offor  EDT, for MDT, and maxfor CDT is determined, where (x', y ) ranges over row R y .2. The First PhaseThe object of the first phase is to determine the function G . We first observethat we can split G into two functions GT (top) and GB (bottom), such thatmin whereWe start with the computation of GT by introducing an array g to store itsvalues. 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 computeg [ x , y ] := GT(x , y) using only g[x , y – 1] in a simple column scan from top tobottom. Similarly, we find GB(i , y ) = GB(i ‚ y + 1) + 1. The second scan runsfrom bottom to top, and computes G(x , y ) directly, using GT from the previousscan, and GB from the current one. After some simplification, this results inthe code fragment given in Fig. 1. Clearly, the time complexity is linear in theA GENERAL ALGORITHM FOR COMPUTING DISTANCE TRANSFORMS 333Fig. 1. Program  fragments  for the first phase.number of pixels (i.e. O (m × n )). In actual implementations it is convenientto replace ∞ by m + n , since all distances in the images are less than m + n ifthe set B is non-empty.3. The Second PhaseIn 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 constantand 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 moregeneral 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 minimizationproblem of Eq. (1). For any i with 0 ≤ i < m , denote by F i the functionon the real interval [0, m – 1]. We call i the index of F i . In the caseof EDT, the graph of F i is a parabola with vertex at (i ‚ g ( i )). In the case of334 A. MEIJSTER ET AL.(a) EDT (b) MDT (c) CDTFig. 2. DT as the lower envelope (solid line) of curves F i , 0 ≤ i < m (dotted lines). Thedashed vertical  lines  indicate  the  transitions  between regions.MDT the parabolas are replaced by V-shaped approximations, while in the caseof CDT we deal with ‘topped off’ V-shaped approximations (see Fig. 2). Wecan interpret DT geometrically as the lower envelope of the collection {F i ⏐ 0 ≤i < m } evaluated at integer coordinates, cf. Fig. 2. The lower envelopesconsist of a number of consecutive curve  segments, whose index we denote bys[0]‚ s [1]‚..., s [q] counting from left to right. The projections of the segmentson the x-axis are called regions , and form a partition of the interval [0, m ) byconsecutive segments. The computation of DT now consists of two scans. In aforward (left-to-right) scan the set of regions is determined using an incrementalalgorithm. In a backward (right-to-left) scan the values DT(x , y ) are triviallycomputed for all x.We start by replacing the upper bound m in (1) by a variable u and defineThe geometric interpretation is that we restrict the set B to the half plane tothe 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 atx 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 hwith 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 scanfrom left to right, and the sets T (h, u ) of points with the same least minimizerh. We thus define(3)A GENERAL ALGORITHM FOR COMPUTING DISTANCE TRANSFORMS 335(a) above (b) below (c) intersectionFig. 3. Location of F u (dashed curve) w.r.t. the lower envelope (solid line).Clearly, S (u  is a nonempty subset of [0, u ), and We)define the regions for u to be the sets T ( h , u ) that are nonempty. It is easy tosee that the regions for u form a partition of [0, m ).The aim is the case where u = m . Indeed, for x ∈ T (h , m ), we haveH (x, m ) = h and hence DT(x , y ) = f (x , h ). The second phase of the algorithmtherefore consists of two scans: scan 3 computes the partition of [0, m ) thatconsists of the regions for m and scan 4 uses these regions to compute DT. Forgiven u, only the curves with indices from 0 to u – 1 are taken into account. Theminimizer of x corresponds to the index of the curve segment whose projectionon the horizontal axis contains x. Let the current lower envelope consist ofq + 1 segments, i.e. S(u) = {s[0], s[1], . . . , s [q]}, with s[l] the index of the l-thsegment. Consider what happens when F u is added. Three situations mayoccur:(a)(b)(c)F u is above the current lower envelope on [0, m – 1], cf. Fig.  ThenS (u + 1) = S (u), since the set T ( u , u + 1) is empty.F u is below the current lower envelope on [0, m – 1], cf. Fig.  ThenS (u + 1) = {u}, i.e., all old regions have disappeared, and there is one newregion T ( u , u + 1) = [0, m ).F u intersects the current lower envelope on [0, m – 1], cf. Fig.  Thecurrent regions will either shrink or disappear, and there is one new regionT (u , u + l).We start searching from right to left for the current region which is inter-sected by F u . This can be determined by comparing the values of F u and F l atthe begin point t[l ] of each current region l = q , q – 1, . . ., until we find the firstl = l * such that Then F u is not the least minimizerat t [ l *], and there must be an intersection of F u with F l * in region l*. Let x *be the horizontal coordinate of the intersection. If l* = q and x * ≥ m we havecase (a); if l* < 0 we have case (b); otherwise case (c) pertains.336 A. MEIJSTER ET AL.Fig. 4. Program fragments for the second phase.To find x *, we introduce a function Sep, where Sep(i , u ) is the first integerlarger or equal than the horizontal coordinate of the intersection point of F uand Fi with i < u , i.e.(4)We thus have x* = Sep(s[ l *], u ). Clearly, the function Sep is dependent onwhich distance transform we want to compute. In the next section we willderive the expressions for the function Sep, but in the remainder of this sectionwe simply assume that Sep is available.We introduce an integer program variable u . It is convenient to representS(u ) by an increasing sequence of elements. Since the regions form a partitionof [0, m ) by consecutive segments, we can represent them by the sequence oftheir least elements. According to the case analysis above, the regions are tobe adapted at their end. We can therefore implement these sequences in twointeger arrays, s and t, with an integer variable q as index of the end point.We start with the forward scan, see scan 3 in Fig. 4. We have S(1) = {0},and T (0, 1) = [0, m ), and thus start with q = 0, s [0] = 0, and t[0] = 0. In aloop, variable u is incremented, and thus the representations of S and T mustbe updated by means of the case analysis above. For details, we refer to ourreport [4].To investigate the complexity of the forward scan, we consider the expressionq + 2(m – u ), which is initially 2m . In every execution of the body of the outerloop (scan 3 in Fig. 4), and also in every execution of the body of its inner loop,A GENERAL ALGORITHM FOR COMPUTING DISTANCE TRANSFORMS 337Fig. 5. Cases for finding  Sep for  MDT.the value of the expression decreases. This implies that the time complexityof the scan is linear in m . Note that, the average number of iterations of theinner loop is at most two. The algorithm uses less than 2m comparisons of ƒvalues, and function Sep is evaluated less than m times.When the forward scan is finished, we have completely determined the parti-tion of [0, m) in regions. Given these regions, we can trivially compute dt-valuesin a simple backward scan (see scan 4 in Fig. 4).4. Derivation of the Function SepThe derivation in the previous section was independent of the actual metricused The functions dependent on the metric are ƒ and Sep. In this section we.compute expressions for Sep for EDT, MDT, and CDT. The easiest is EDT.We find for i < uHere, we denote integer division with rounding off towards zero by div. Thus,we find for EDT thatIf we use the Manhattan metric, the analysis is slightly more complicated.Since we have to deal with absolute values in the expressions, awkward caseanalysis is necessary if we want to compute Sep analytically. Therefore weprefer a geometric argument. We have to consider three cases (see Fig. 5).If g(u) ≥  g ( i ) + u – i , the graph of F u lies entirely above the graph of Fifor all x, thus we choose Sep( i, u ) = ∞ . If g (i) > g (u) + u – i , the graph of F ilies entirely above the graph of F u , so Fi (x ) ≤ F u (x ) for no x at all. Thus, wemust choose Sep( i, u ) = –∞  to satisfy (4). In all other cases, F u intersects F i338 A. MEIJSTER ET AL.Fig. 6. Cases for finding Sep for CDT, where γ = (u – i )/2. CasesCasesat So, if we want to compute MDT we useFor the case of CDT we have |x – i | max g (i ) ≤ |x – u | max g (u). Weconsider two main cases, which each can be split up in three sub-cases. Firstwe consider the case g(i) ≤ g (u). From Fig.  we see that the increasingsegment of F i (y = x – i) intersects the decreasing part of F u ( y = u – x ), orthe constant part (y = g (u)). Let γ be the vertical coordinate correspondingwith the middle of i and u (x = (i + u )/2), i.e. γ = (u – i )/2. From Fig. we see that if we haveFrom Fig. 6(b)-(c), we see that the increasing part of F i intersects the constantsegment of F u at i + g (u), and thus we havePutting the three cases together, we can concludeA GENERAL ALGORITHM FOR COMPUTING DISTANCE TRANSFORMS 339TABLE ITiming results in ms. From left to right: EDT, MDT, and CDT.The other main case is g(i) > g (u). Again, in Fig.  we see that if g (i) ≤ γ,the intersection at (i + u )/2 is the separator. If g( i) > γ (see Fig.  thehorizontal segment of F i intersects the decreasing part of Fu at x = u – g ( i).Just like in the previous case, we can put these cases together. This results inthe following expression for Sep:5. Parallelization, Timing Results, and ConclusionsSince the computation per row (column) is independent of the computationof other rows (columns), the algorithm is well suited for parallelization on ashared memory machine. In the first (second) phase, the columns (rows) aredistributed over the processors. The two phases must be separated by a barrier ,which assures that all processors have completed the first phase before any ofthem starts with the second phase. The theoretical time complexity of theparallel algorithm for p processors (where p ≤ m min n ) is (m n/ p).We ran experiments on an Intel Pentium III based shared memory parallelcomputer with 4 cpu's, running at a 550MHz clock frequency. We performedtime measurements using several binary images, and found that the executiontime is almost independent of image content, and scales well w.r.t. the numberof processors. This is as expected, since the amount of work per row andcolumn is almost the same. In table I the timings for square images are givenfor p = 1 to p = 4 processors. Note that the computation of MDT and CDT isonly slightly faster than the exact EDT. We also implemented the sequentialalgorithm of [7] for CDT, and found that our algorithm is less than a factor of2 slower, which can easily be overcome by parallel processing.The algorithm can be easily extended to d -dimensional distance transformsby separating the problem into d phases, each solving a one-dimensional prob-lem, as carried out above for the case d = 2.340 A. MEIJSTER ET AL.References1 .2 .3 .4 .5 .6 .7 .G. Borgefors. Distance transformations in arbitrary dimensions. Computer Vision,Graphics, and Image Processing, 27:321–345, 1984.G. Borgefors. Distance transformations in digital images. Computer Vision, Graphics,and Image Processing, 34:344–371, 1986.P. Danielsson. Euclidean distance mapping. Comput. Graphics Image Process., 14:227–248, 1980.W. H. Hesselink, A. Meijster, and J. B. T. M. Roerdink. An exact Euclidean distancetransform in linear time. Technical Report IWI 99-9-04, Institute for Mathematics andComputing Science, University of Groningen, the Netherlands, Apr. 1999.M. Kolountzakis and K. Kutulakos. Fast computation of the Euclidean distance mapsfor binary images. Information Processing Letters, 43:181–184, 1992.S. Pavel and A. Akl. Efficient algorithms for the Euclidean distance transform. ParallelProcessing  Letters, 5:205–212, 1995.A. Rosenfeld and J. Pfaltz. Distance functions on digital pictures. Pattern Recognition,1:33–61,  1968.

A general algorithm for computing distance transforms in linear time

https://pure.rug.nl/ws/files/3059926/2002CompImagVisMeijster.pdf

A general algorithm for computing distance transforms in linear time

Abstract

Similar works

Full text

Available Versions

ARTS repository - University of Groningen