Fast computation of the Euclidian distance maps for binary images

Information Processing Letters 43 (1992) 181-184 North-Holland

28 September 1992

Fast computation of the Euclidian distance maps for binary images Mihail N. Kolountzakis

*

Department of Mathematics, Stanford Unit:ersity, Stanford, CA 94305, USA

Kiriakos N. Kutulakos

*

Computer Sciences Department, Unit,ersity of Wisconsin - Madison, Madison, WI 53706, USA Communicated by D. Dolev Received 15 December 1991 Revised 13 July 1992

Abstract Kolountzakis, M.N. and K.N. Kutulakos, Fast computation of the Euclidian distance maps for binary images, Information Processing L e n e r s 43 (1992) 181-184. A simple algorithm is given for the computation of the Euclidian distance from the set of black points in an N X N black and white image, for all points in the image. The running time is O ( N 2 log N ) and O ( N ) extra space is required. The algorithm is suitable for implementation on a parallel machine.

Keywords: Euclidian distance map; parallel algorithms; computer graphics; pattern recognition; robotics; image processing

1. Introduction Consider a black and white N X N binary image: i.e., a two-dimensional array where aij = 0 or 1, for i,j = 0 , . . . , N - 1 . The index i stands for the row, the index j for the column and (0, 0) is the upper left point in the image. We solve the problem of finding for each point (i, j) its Euclidean distance from the set of all black points B = {(i, j): a u = 1}. In other words, we compute the array

d~=

rain

((i-x)2+(j-y)2),

foralli,j.

(x,y)cB -

Correspondence to: M.N. Kolountzakis, Department of Mathematics, Stanford University, Stanford, CS 94305, USA. * Both authors were at the Image Analysis Laboratory, Department of Computer Science, University of Crete, Greece, when this work was carried out (1988).

An algorithm is given which requires O ( N 2" log N ) time and O ( N ) space (plus the space needed to store the result). The algorithm is suitable for implementation on a parallel machine. The Euclidean distance map d u has important uses in computer vision, pattern recognition and robotics [2]. For instance, if the black points represent obstacles then d u tells us how far the point (i, j) is from these obstacles. This information is useful when one tries to move a robot in the free space (white points of the image) keeping it away from the obstacles. It is much easier to compute the distance map for the L ~ distance. This can be done in optimal time O ( N 2 ) , if one uses the fact that the L ~ distance between two points in the N × N grid is the length of the shortest path that joins them. However, there is no such interpretation for the

0020-0190/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved

181

Volume 43, Number 4

I N F O R M A T I O N PROCESSING LETTERS

Euclidean distance. O n e way a r o u n d this problem is to a p p r o x i m a t e the Euclidean distance with distances that can be c o m p u t e d in the same way as the L l distance [2]. T h e only nontrivial algorithm for the c o m p u t a t i o n of the exact Euclidean distance m a p that we are aware of is that of Y a m a d a [3]. Y a m a d a ' s algorithm p e r f o r m s N iterations of local t r a n s f o r m a t i o n s of the entire image. It is not clear w h e t h e r Y a m a d a ' s algorithm can be converted to an efficient (i.e. running in time close to O ( N 2 ) ) algorithm on a sequential c o m p u t e r .

2. The algorithm First note that it is sufficient to give an algorithm for the c o m p u t a t i o n of r~ =

rain ((i-x) (x,y)~B,y>~j for all i,j.

2 + (j-y)2),

H e r e r 2. is the square of the distance of the point (i, j) from the part of B that is to the right of (i, j'). C o m p u t i n g the distance from the right and from the left and taking the m i n i m u m of the two gives d/Zj. So only the c o m p u t a t i o n of rij is described here. Let Bj = B C3 {(x, y): y > j } and let N,j denote any n e a r e s t point to (i, j) in Bj. T h e following l e m m a is exploited by our algorithm.

Lemma 1. Let P = ( a , j ) , Q = ( b , j ) , b < a , be two points in the same column with Q above P. Let Np = ( x, y ), and N o = ( z, w ). Then z 4 x, that is, NO is above Np. This holds for any choice of the nearest black points N e and NQ. Proof. W e have (x-

a) 2 + ( y _ j ) 2 4 ( z _ a )

2 + (w-j)2

and (z-b)2+(w-j)2~

(x-b)

~+ (y-j)2.

Adding the above inequalities gives z ~
[]

T h e array N~i is c o m p u t e d column by column, moving from left to right. T h e values of r~2 can be 182

28 September 1992

c o m p u t e d from N~j. While computing fixed j and for all i, we use the vector

Nij for

V/°) = m i n { k : k > j , (i, k ) e B } . If the above set is empty, we set Vi°) = o0. Each vector Vii°+~) is c o m p u t e d from Vi°) as follows. For each i such that aij = 1, scan the ith row to the right, starting from column j + 1, until a black point is e n c o u n t e r e d , which will be V/i°+1). If no black point is e n c o u n t e r e d up to the column N - 1, set V/(j+l) = o0. If air = 0 then set Vi(j+~) = Vi (j).

Only the current V °) vector is stored. T h e c o m p u t a t i o n for all i and j of Vi°) takes O ( N 2) time, since each row is traversed to the right only once. So, while c o m p u t i n g N~j for any fixed colu m n j, we may assume the vector V °) is given. It remains to give an O ( N log N ) time algorithm for c o m p u t i n g Nit for any fixed j. We have to find for each i a nearest point of (i, j) a m o n g the points in the vector V °). (Notice that we identify V~j) with the point (k, Vk~J)).) Doing this by exhaustive search results in an O ( N 2) algorithm. Instead, the following recursive algorithm is used. D e n o t e by A ( x 1, x2, zl, z 2) the algorithm below. This algorithm finds nearest black points for all points (i, j), x I ~< i ~
Algorithm A(xl, Xz, Zl, z 2) 1

1. Set x 0 = [7(x I +x2)]. 2. Find N~d = (z0, w 0) by searching one by one the points V2 i), z 1 ~< k ~< z 2. 3. If Xl ~
Volume 43, Number 4

INFORMATION PROCESSING LETTERS

28 September 1992 4

R e m a r k on other distances

where also Y'. "i ,~(2)_ -n+1+2

O u r algorithm can be applied to the c o m p u t a tion of the distance map for distances other than the Euclidean distance. Such a distance must satisfy L e m m a 1 and must also have the following property: if b ~
i=1

L

_
1 2k

(m

2L

E k=l

i=__

i=1

n?'

)

2L '

where for all k, 2k E ~ ("i k)

=n+l+2+

'"

+2k-~
i=1

2 l-

~
T ( I , n(i")).

i=1

H e n c e there is a constant C > 0 such that

T ( r n , n) <~ C ( n + m + (rn + n ) log m ) . This gives the O ( N log N).

desired

T(N-1,

N-1)=

4. Parallelizing the algorithm 3. Time Let T ( m , n) be the maximum running time of algorithm A when x 2 - x ~ <~m and z 2 - z ~ ~ n . Obviously, T(1, n ) = n and in general we have the recurrence inequality

T(m,n)<~max

n+T

,n I +T

,n 2

,

where n~ + n 2 = n + 1, and the maximum is taken over all such nl,n 2. The first two steps of the algorithm account for the n term and the two recursive calls account for the two T terms in the right-hand side. A s s u m e for simplicity that m = 2 c is a power of 2. Iterating the above inequality L times we have

T ( m , n) <~n+T

, n I) + T

,n

where n] l) + n~21)= n + 1 ~n

-~- n ~ l ' - [ - n(21) q- £ i=1

r

, n{ 2)

Since the above algorithm works on one column at a time, it can be easily parallelized. Consider for example an Exclusive Read, Exclusive Write Parallel R a n d o m Access Machine ( E R E W P R A M ) with p processors [1]. This consists of r a n d o m access m e m o r y to which all p processors have simultaneous access provided they read or write to different locations. O u r algorithm can be parallelized on the E R E W P R A M with optimal speedup as follows. We use one more N × N array to store the values of Vi(j) for all i and j. To c o m p u t e all V/i~j) we assign N / p conscecutivc rows to each processor. Each processor is responsible for computing the values of Viii) for all i and j that belong to the region that has been assigned to the processor. This does not give rise to any r e a d / w r i t e conflicts since all processors work on disjoint regions. This c o m p u t a t i o n takes O(NZ/p) time, since each row of the V array takes O ( N ) time, and each processor c o m p u t e s N / p rows. Having c o m p u t e d and stored the values of ~(i), we now assign N / p consecutive columns to each processor. Each processor can perform the algorithm A(0, N - 1, 0, N - 1) on all the 183

VoLume 43, Number 4

INFORMATION PROCESSING LETTERS

c o l u m n s that have b e e n a s s i g n e d to it, since, for each column, all the i n f o r m a t i o n that the p r o c e s sor n e e d s is s t o r e d in the c o r r e s p o n d i n g c o l u m n of the V array. A g a i n no r e a d / w r i t e conflicts arise, the r u n n i n g time is O ( N 2 log N / p ) a n d the s p e e d u p is p. A s m e n t i o n e d in Section 2 o u r a l g o r i t h m applies w i t h o u t any c h a n g e s to the case of the L 1 distance. O f course, for the L ~ d i s t a n c e t h e r e is a well known a n d o p t i m a l O ( N 2) algorithm. This a l g o r i t h m starts p r o p a g a t i n g a wave from the black points until it has r e a c h e d the w h o l e image. This a l g o r i t h m is s u p e r i o r to ours in the case of the s e q u e n t i a l i m p l e m e n t a t i o n , but a p a r a l l e l imp l e m e n t a t i o n of it is far f r o m obvious, since the access to the i m a g e is not uniform. In this case o u r a l g o r i t h m could be p r e f e r a b l e .

184

28 September 1992

Acknowledgment W e w o u l d like to t h a n k the r e f e r e e s for their helpful c o m m e n t s .

References [1] T.H. Cormen, C.E. Leiserson and R.L. Rivest, Introduction to Algorithms (McGraw-Hill, New York, 1990). [2] P.E. Danielsson, Euclidean distance mapping, C'omput. Graph. hnage Processing 14 (1980) 227-248. [3] H. Yamada, Complete Euclidean distance transformation by parallel operation, in: Proc. 7th Internat. ('onj~ o;; Pattern Recognition. Vol. I (1984) 69 71.