On Digital Distance Transforms in Three Dimensions

COMPUTER VISION AND IMAGE UNDERSTANDING

Vol. 64, No. 3, November, pp. 368–376, 1996 ARTICLE NO. 0065

On Digital Distance Transforms in Three Dimensions GUNILLA BORGEFORS* Centre for Image Analysis, Swedish University of Agricultural Sciences, La¨gerhyddva¨gen 17, 752 37 Uppsala, Sweden Received November 30, 1994; accepted August 2, 1995

Digital distance transforms in 3D have been considered for more than 10 years. However, not all of the complexities involved have been unravelled. In this paper the complete geometry and equations for 3D transforms based on a 3 3 3 3 3 neighborhood of local distances are given. A new type of valid distance transforms (DTs) have been discovered. The optimal solutions are computed, where optimality is defined as minimizing the maximum difference from the true Euclidean distance, thus making the DTs as direction independent as possible. The well-known k3, 4, 5l DT is confirmed as the most practical weighted DT, where the distance is set to 3 between neighbors sharing an area, 4 between neighbors sharing an edge, and 5 between neighbors sharing a point.  1996 Academic Press, Inc.

1. INTRODUCTION

Digital shape is an important and much used concept in image analysis. One excellent tool for many different operations related to shape are digital distance transforms. In a distance transform (denoted DT), each element in the shape (background) has a value measuring the distance to the background (shape). DTs can, e.g., be used for morphological operations, skeletonization, computing Voronoi diagrams, and template matching. For an overview of applications, see [6]. The basic idea, utilized for most DTs, is to approximate the global Euclidean distance by propagation of local distances, i.e., distances between neighboring pixels. This idea was probably first presented by Rosenfeld and Pfaltz in 1966 [10]. This approach is motivated by ease of computation. In sequential computation only a small area of the image is available at the same time. In massively parallel computation each pixel has access only to its immediate neighbors. Now that 3D voxel images are common, not only in medical imaging but also in other areas, it is useful to consider distance transforms in 3D digital space. In fact, the first results in this area was published over 10 years ago, and there have been a number of papers since then. * E-mail: [email protected]. 368 1077-3142/96 $18.00 Copyright  1996 by Academic Press, Inc. All rights of reproduction in any form reserved.

Danielsson, in 1980, published the first paper on how to reasonably compute the true Euclidean distance transform in digital images [7]. He included the 3D case. Ragnemalm has recently improved these algorithms [9]. Okabe, Toriwaki, and Fukumura in 1983 extended the city block/chessboard and octagonal DTs to 3D [8]. In these DTs, only the number of steps necessary to connect one voxel with another is considered, where connectivity is dependent on the definition of neighborhood. Borgefors in 1984 considered many types of DTs in 3D, including weighted ones [3]. The maximum difference from the Euclidean distance in an M 3 M 3 M image was computed and minimized by finding the optimal local step lengths. A slightly different approach is that by Verwer, 1991 [12]. He minimizes the maximum difference, not in a cubic M 3 M 3 M image, but on a (nondigital) Euclidean sphere with radius M, embedded in digital space. A related problem, that has also given insights in this area, is the approximation of straight lines in digital space. Beckers and Smeulders, in 1992, minimized the mean squared difference between the computed line length and the Euclidean line length [2]. The optimal values they find for the local step lengths can be used with good results as weights when computing DTs. Note, however, that the problem of finding the best digital approximation for a random straight line in real space is not the same as finding the best DT. In the latter case there are no underlying nondigital considerations. The three last above-mentioned works [3, 12, 2] also consider integer-valued approximations of the optimal local step lengths. With all this effort gone into 3D distance transforms, one would think that there is very little to add on the subject. However, 3D space is more topologically complex than previously realized. It is often impossible to simply extend 2D results. In the present case the difficulty arises both from the fact there are so very many paths connecting one voxel to another in digital space and from the fact that the solutions of the optimization procedures used are valid in only very limited regions of parameter space. Only Beckers and Smeulders [2] have so far considered the limits of their solutions at all. In this paper the geometry, local distance value limits,

369

DIGITAL DISTANCE TRANSFORMS

and equations of 3D distance transforms based on a 3 3 3 3 3 voxel neighborhood will be completely described. Previous efforts (including those by the present author) have found only one type of DTs and missed the second, equally valid, type. The results are not only of value for computing DTs, but should contribute to the knowledge of 3D digital space in general. In Section 2 the geometry and general equations are developed. This is the central part of the paper. In Section 3 optimal real and integer-valued DTs are computed, where optimality is defined as minimizing the maximum difference from the Euclidean distance in an M 3 M 3 M image. Readers preferring other optimality criteria, or not interested in lots of mathematical formulae, can skip this section. In Section 4 the optimal weighted DTs are listed and the two basic types of DTs are illustrated by their associated spheres. Finally, Section 5 briefly reminds the reader how to compute 3D DTs. 2. GEOMETRY AND EQUATIONS

Denote a digital shape on a cubic grid F, and the complement of the shape F , where the sets F and F are not necessarily connected. A distance transformation converts the binary image to a distance image, or distance transform. In the DT each voxel has a value measuring the distance to the closest pixel in F . A good underlying concept for all digital distances is the one below proposed by Yamashita and Ibaraki, [14]. DEFINITION 1. The distance between two points x and y is the length of the shortest path connecting x to y in an appropriate graph. They proved that any distance is definable in the above manner, by choosing an appropriate neighborhood relation and an appropriate definition of path length. In 3D cubic space each voxel has three types of neighbors: six area neighbors, 12 edge neighbors, and eight point neighbors. A path between two voxels in the 3D image can thus include steps in 26 directions, if only steps between immediate neighbors are allowed. The distance transform is defined from the minimal paths as follows. DEFINITION 2. The distance transform DT(i, j, k) of a binary image is defined as DT(i, j, k) 5 min hD(maxhui 2 pu, u j 2 qu, uk 2 ruj, (p,q,r)[F

midhui 2 pu, u j 2 qu, uk 2 ruj, minhui 2 pu, u j 2 qu, uk 2 ruj)j,

where F is the set from which the distances are measured, mid is the middle of the three values, and D(?, ?, ?) is the distance between two pixels defined by a minimal path between them.

Due to symmetry, it is enough to consider distances from the origin to a voxel (x, y, z), where 0 # z # y # x # M and M is the maximal dimension of the image. The distance to be minimized in Definition 2 then simply becomes D(x, y, z). In path generated distance transforms (PGDT) the distance is defined as the number of steps in the minimal path. The three simplest 3D path-generated distances are the ones where steps can be taken to the six closest neighbors (area neighbors); to the 18 closest neighbors (area 1 edge neighbors); and to the 26 closest neighbors (area 1 edge 1 point neighbors). These DTs are denoted D6, D18, and D26. D6 can be considered the 3D equivalent of the city block distance (D4) and D26 as the 3D equivalent of the chessboard distance (D8). Just as the city block and chessboard distances can be mixed in 2D to create octagonal DTs that are much better approximations to the Euclidean distance (see [11]), so can the three basic PGDTs in 3D be mixed for better approximations. These issues were investigated by Okabe, Toriwaki, and Fukumura [8] and independently discovered and illustrated in [3]. The disadvantage with these ‘‘mixed’’ solutions is that they are computationally intensive, unless a massively parallel architecture (one processing element per voxel) is available. As the length of any minimal path is defined only by the numbers of steps of different types in it, the order of the steps is arbitrary. Therefore, we can always assume a minimal path where the steps are arranged in a number of straight line segments, equal to the number of different directions of steps used. In Fig. 1 the number of line segments in the minimal path is shown for the three PGDTs. The figure illustrates the minimal paths from the origin to a quarter slice of space at height z 5 5. All voxels with distance value ,10 are marked. (When different minimal paths using different step combinations exist, the one with the least number of directions is chosen.) In D6 only the voxel (0, 0, 5) can be reached by one straight line segment, while in D26 (0, 0, 5), (0, 5, 5), (5, 0, 5), and (5, 5, 5) can all be reached in this way. Contrary to expectations, D18 uses four different directions in many cases, resulting from edge neighbor steps in two different directions. This distance is irregular in many ways, which may explain why several previous authors have wisely mentioned only D6 and D26 [12, 2]. Mathematical expressions for the distance between two voxels are necessary for computing various DT properties. The three PGDTs become (from [8]) D6 5 x 1 y 1 z,

(1)

HK

D18 5 x 1 max 0, D26 5 x,

z1y2x11 2

HJ,

(2) (3)

370

GUNILLA BORGEFORS

FIG. 1. The number of directions in the minimal paths. One quarter slice of space at z 5 5, where the origin is below the top left voxel.

where x, y, z denote the ordered differences between voxel coordinates, as in Definition 2. The PGDTs only count the number of steps, which means that all steps are considered having length one. These DTs are closely related to connectedness. As the local Euclidean distances between the three types of neighbors are 1, Ï2·, and Ï3·, respectively, this leads to very rough approximations of the global Euclidean distance. By using different local distances, the weighted (or chamfer) distance transforms (WDT) are defined. Denote the local length between area neighbors a, between edge neighbors b, and between point neighbors c. A 3 3 3 3 3 WDT will be denoted ka, b, cl. Not all combinations of local distances a, b, c result in useful distance transforms. The following property is desirable: DEFINITION 3. Consider two pixels that can be connected by a straight line, i.e., by using only one type and direction of local step. If that line defines the distance between the pixels, i.e., is a minimal path, then the resulting DT is semi-regular. If there are no other minimal paths, then the DT is regular. In 2D a regular distance transform is a metric, as all regular DTs are positive linear combinations of the city block and the chessboard DTs, which have been proved to be metrics [10, 11]. The situation has not been resolved in 3D, but the conjecture is that regular DTs are metrics here too. As there are three types of steps, there are only three types of straight paths possible in the cubic grid. To find the conditions for 3D regularity we must investigate all the ways these three straight paths can be approximated by paths using other steps and find the conditions for the straight path to be shortest. The result of this investigation is found in Table 1. For example, if a c-line is approximated by b-steps, then three such steps are necessary to replace two c-steps: (1, 1, 0) 1 (1, 0, 1) 1 (0, 1, 1) 5 2 3 (1, 1, 1). For regularity, 2c , 3b must be valid. This is the next to last line in Table 1. Not all of the nine inequalities are significant. The easiest way to see which conditions are necessary and sufficient is to draw the inequalities in a diagram, using b and c as axes and a as a scaling factor;

see Fig. 2. Four inequalities define an area of regularity in parameter space. They are a , b , 2a and b , c , Dsb.

(4)

The conditions in (4) may seem restrictive, but they are not sufficient to determine unique expressions similar to (1)–(3) for the WDTs. If we compute the distances from the origin, choosing the shortest paths and assuming that the local distances have the properties in (4), we discover two different, equally valid cases. They are exemplified at voxel (2, 1, 1). It can be reached either as (1, 1, 1) 1 (1, 0, 0) 5 a 1 c or as (1, 1, 0) 1 (1, 0, 1) 5 2b. The inequalities in (4) do not determine which is the shorter path. Thus, in addition to the regularity conditions in Eq. (4), we have Case I: a 1 c , 2b, Case II: a 1 c . 2b.

(5)

The line a 1 c 5 b dividing the two cases is drawn dashed in Fig. 2. The computed DT values at all voxels 0 # z # y # x # 5 in both cases are found in Figs. 3 and 4. From these figures and some generalization, the weighted distance between the origin and (x, y, z) becomes Case I: D 5 z(c 2 b) 1 y(b 2 a) 1 xa,

TABLE 1 Regularity Criteria for 3 3 3 3 3 Distance Transforms Path

Steps

Inequality

a a a b b b c c c

b c b c a c a c a b a b

2a , 2b 2a , 2c a,b1c b , 2a 2b , 2c b,c1a c , 3a 2c , 3b c,a1b

(6)

DIGITAL DISTANCE TRANSFORMS

371

FIG. 2. The cone in a, b, c-space that results in regular 3D weighted distances. The dashed line separates the two cases. FIG. 4. The distance values from the origin for 3 3 3 3 3 weighted distance transforms Case II for 0 # z # y # x # 5.

Case II: D 5

FIG. 3. The distance values from the origin for 3 3 3 3 3 weighted distance transforms Case I for 0 # z # y # x # 5.

5

z(c 2 b) 1 y(c 2 b) 1 x(2b 2 c) if x # y 1 z

z(b 2 a) 1 y(b 2 a) 1 xa if x $ y 1 z.

(7)

Case I is the only one found in the literature [3, 12, 2]. Beckers and Smeulders [2], who seem to be the only authors to have seriously considered admissible relations between local distances, state without proof that if the restrictions b # 2a and a 1 c # 2b are valid, then the resulting WDT is a metric. If they have also made the ‘‘natural’’ assumption a # b # c these conditions correspond to Case I, as seen in Fig. 2. To further emphasize the difference between the two cases, the number of directions in the minimal paths have been computed, just as in Fig. 1. The results are found in Fig. 5. The number of directions are, of course, the same for all Case I and Case II WDTs, respectively. Here good integer approximations are used to represent the two cases (see Sections 3.3 and 4 for discussions on integer WDTs). Note that in Case II, even though only b-steps may be used in some paths, the path may still have two directions. An example is, of course, (2, 1, 1) 5 (1, 1, 0) 1 (1, 0, 1). Note, also, that no minimal path uses more than three directions.

372

GUNILLA BORGEFORS

and fmax 5 (b 2 Ï1 2 a 2 )M for h 5

aM Ï1 2 a 2

if j ; 0.

Proof. The extremal value is found by setting the derivative of f 9(h) to zero (regarding all other parameters as constants), solving for h, and simplifying the resulting expressions. n FIG. 5. The number of directions in the minimal paths for Case I (left) and Case II (right) weighted distances. Integer WDTs of both types are used as examples. One quarter slice of space at z 5 5, where the origin is below the top left voxel (cf. Fig. 1).

It is possible to consider the PGDTs D6 and D26 as weighted DTs, by choosing the values of the three steplengths appropriately. We get D6 p k1, 2, 3l and D26 5 k1, 1, 1l. Both DTs are semiregular. For D18 the situation more complex. From Eq. (2) we have

cp

5

3 2

if z even,

3 1 1 if z odd. 2 2z

3. OPTIMALITY COMPUTATIONS

In this section the optimal local distances for the two regular 3D WDT cases will be computed. Optimality is defined as minimizing the maximal difference between the WDT and the Euclidean distance in a cubic image of size M 3 M 3 M. The choice of optimality criteria is somewhat arbitrary. However, this one has the advantage that it does not depend on any nondigital structure, such as an embedded Euclidean sphere. The maximum of a particular type of function will often have to be computed. The following lemma is used.

f (h) 5 ah 1 cj 1 bM 2 ÏM 2 1 j 2 1 h2, where uau , 1 and M . 1 has the maximum value fmax 5 bM 1 cj 2 ÏM 2 1 j 2 ? Ï1 2 a 2 for h 5

aÏM 2 1 j 2 if j ? 0 Ï1 2 a 2

The difference between the computed distance (see (6)) and the Euclidean distance is Diff I 5 (c 2 b)z 1 (b 2 a)y 1 ax 2 Ïx 2 1 y 2 1 z 2,

(9)

where 0 # z # y # x # M. This difference is to be minimized. Let x 5 M. Then 0 # y # M and 0 # z # y. The extrema for Diff I(z) occur for z 5 0; ­/­z(Diff I) 5 0; or z 5 y. Using Lemma 1 for the middle value, while assuming h 5 z, j 5 y, a 5 c 2 b , 1, c 5 b 2 a, and b 5 a, and simple insertion in (9) for the other two, the possible maxima become

(8)

Thus, D18 p k 1, 1, Ds 1 «l, which is not regular, as « $ 0.

LEMMA 1. The function

3.1. Case I

E I1,y 5 (b 2 a)y 1 aM 2 ÏM 2 1 y 2 E I2,y

for z 5 0,

5 (b 2 a)y 1 aM 2 Ï1 2 (c 2 b) ÏM 1 y

E I3,y 5 (c 2 a)y 1 aM 2 ÏM 2 1 2y 2

2

2

2

for 0 , z , y, for z 5 y.

For each of these three error expressions the maximum can occur for y 5 0; ­/­y(E Ii ) 5 0; or y 5 M. Simple insertion and Lemma 1 (with j 5 0) give seven possibilities, if we assume b 2 a , 1. Note that E I1,1 5 E I2,1 5 E I3,1 , as z # y: E I1,1 5 (a 2 1)M

for (M, 0, 0),

E I1,2 5 (a 2 Ï1 2 (b 2 a)2 )M

50 for (M, y zmax , 0),

E I1,3 5 (b 2 Ï2)M

for (M, M, 0),

,y ,y E I2,2 5 (a 2 Ï1 2 (c 2 b)2 2 (b 2 a)2 )M for (M, y zmax , z zmax ),

E I2,3 5 (b 2 Ï2 Ï1 2 (c 2 b)2 )M

5M for (M, M, z ymax ),

E I3,2 5 (a 2 Ï1 2 sA (c 2 a)2 )M

5y 5y for (M, y zmax , y zmax ),

E I3,3 5 (c 2 Ï3)M

for (M, M, M). (10)

It can be shown with moderate effort that when c 1 z,y z,y a , 2b, which is true in Case I, then y max , z max in E I2,2 . In the other six cases it is obvious that z # y. Numerical experimentation shows that maxhE Ii, j j is minimal for

DIGITAL DISTANCE TRANSFORMS

373

* 5 1, a Iopt b Iopt * 5 Ah (2Ï3 1 2 1 Ï8 Ï Ï3 2 1) P 1.31402, (13) * 5 Ad (2Ï3 2 1 1 Ï8 Ï Ï3 2 1) P 1.62803, c Iopt maxdiff I* 5 Ad (Ï3 1 1 2 Ï8 Ï Ï3 2 1) P 0.10402. (14) 3.2. Case II The difference between the computed distance (see (7)) and the Euclidean distance is Diff eII 5 (z 1 y)(c 2 b) 1 x(2b 2 c) 2 Ïx 2 1 y 2 1 z 2 if x # y 1 z, Diff kII 5 (z 1 y)(b 2 a) 1 xa 2 Ïx 2 1 y 2 1 z 2 if x $ y 1 z. FIG. 6. The difference from the Euclidean distance for the optimal local distances at z 5 100, every tenth voxel included. The maximal positive and negative error are marked in bold italic; the theoretical absolute maximum is 7.36.

(15)

E I2,2 5 2E I1,1 5 2E I1,3 5 2E I3,3 . Solving these equations yields

Again, let x 5 M. Then 0 # y # M and 0 # z # y. We must separate z 1 y # M and z 1 y $ M. The interesting areas can be found in Fig. 7. The maximum of Diff II can occur within the two triangular areas, on the boundaries, or at the corners. Starting with the two areas, k2 and e2 , the extrema occur for ­/­z(Diff II) 5 0. Using Lemma 1 again we get

a Iopt 5 As (1 1 Ï2Ï2(1 1 Ï3) 2 7)

2 2 2 Diff II k2,max 5 y(b 2 a) 1 Ma 2 Ï1 2 (b 2 a) ÏM 1 y ,

P 0.92644,

0 , z , M 2 y, 0 , y , M,

b Iopt 5 As (2Ï2 2 1 1 Ï2Ï2(1 1 Ï3) 2 7) P 1.34065, (11) c Iopt 5 As (2Ï3 2 1 1 Ï2Ï2(1 1 Ï3) 2 7) P 1.65849, with maxdiff I 5 As (1 2 Ï2Ï2(1 1 Ï3) 2 7) P 0.07356. (12) The above optimization only takes the maximum difference between the computed distance and the Euclidean distance into account. The distribution of differences in different directions is also interesting. In Fig. 6 the differences have been computed in a quarter slice of space at z 5 100. Every tenth voxel is listed. The maximum negative difference (equivalent to E I1,1 , E I1,3 , and E I3,3 in (10)) is found in the outer corners and the maximum positive difference (equivalent to E I2,2 in (10)) is found almost in the middle of the slice. The positive differences dominate, which indicates that minimizing the mean difference would not give the same optimal local distances. Figure 6 is, of course, also a check that the above computations are correct. The optimal solutions for a ; 1 will be needed later. In this case we solve E I2,2 * 5 E I3,2 * 5 2E I3,3 * , where the star denotes that a 5 1 has been subsituted in the expressions. The solutions are

Diff e2,max 5 y(c 2 b) 1 M(2b 2 c) 2 Ï1 2 (c 2 b)2 ÏM 2 1 y 2, II

M 2 y , z , y, M/2 , y , M,

where we also assume c 2 b , 1/Ï2· and b 2 a , 1/Ï2·. On the lines in Fig. 7 we have the following expressions for Diff II:

FIG. 7. The different areas, lines, and points where the maximal differences from the Euclidean distance can occur in Case II (see the text).

374

GUNILLA BORGEFORS

Diff kII1 5 (b 2 a)y 1 aM 2 ÏM 2 1 y 2,

0 , y , M,

Diff kII3 5 2(b 2 a)y 1 aM 2 ÏM 2 1 2y 2,

0 , y , M/2,

Diff lII2 5 bM 2 Ï2 ÏM 2 2 My 1 y 2,

M/2 , y , M,

Diff eII1 5 (c 2 b)z 1 bM 2 Ï2M 2 1 z 2,

0 , z , M,

Diff eII3 5 2(c 2 b)y 1 (2b 2 c)M 2 ÏM 2 1 2y 2, M/2 , y , M.

The maxima of all these expressions can be found using Lemma 1 with j 5 0. Together with the values at the corner points in Fig. 7 the possibilities for the maximal difference in Case II are

solution is a worse approximation to the Euclidean distance than the one from Case I (maxdiff I , maxdiff II). It can also be remarked that the value of E lII3 is very, very close to maxdiff II. As in Case I we need the optimal solutions for a ; 1. In this case maxdiff depends only on b, determined from E kII3* 5 2E lII1*, while the local distance c can vary within an interval. The minimum c value is determined by E eII4* and the maximum by E eII3*. The solutions become II *51 a opt II * 5 Ad (Ï2 1 1 1 Ï8Ï2 2 9) bopt

P 1.31178, (19)

E kII1 5 (a 2 Ï1 2 (b 2 a)2 )M

1 for (M, y kmax , 0),

II c opt,min * 5 Ad (3Ï3 2 2Ï2 1 1 1 Ï8Ï2 2 9) P 1.62960,

E kII2 5 (a 2 Ï1 2 2(b 2 a)2 )M

2 2 ), for (M, y kmax , z kmax

II c opt,max * 5 Al (13 2 2Ï2 1 Ï8Ï2 2 9)

E kII3 5 (a 2 Ï1 2 2(b 2 a)2 )M

3 3 ) for (M, y kmax , y kmax

E kII4 5 (a 2 1)M

for (M, 0, 0),

E lII1 5 (b 2 Ï2)M

for (M, M, 0),

E lII2 5 B

as Diff lII2 is monotonous,

E lII3 5 (b 2 ÏsD )M

for (M, M/2, M/2),

E eII1 5 (b 2 Ï2 Ï1 2 (c 2 b)2 )M

1 ), for (M, M, z emax

maxdiffII* 5 Ad (2Ï2 2 1 2 Ï8Ï2 2 9)

3 , y e3 , 0), E eII3 5 (2b 2 c 2 Ï1 2 2(c 2 b)2 )M for (M, y emax max

for (M, M, M). (16)

2 2 3 2 2 3 , thereIn fact, z kmax 5 y kmax 5 y kmax and z emax 5 y emax 5 y emax II II II II fore E k2 5 E k3 and E e2 5 E e3 , respectively. Thus, there are eight expressions to minimize, rather than 11. Numerical experimentation shows that maxhE nIIj j is minimal for E kII3 5 2E kII4 5 2E lII1 5 E eII3 . The solutions are

II 5 As (1 1 Ï4Ï2 2 5) a opt

P 0.90523,

II 5 As (2Ï2 2 1 1 Ï4Ï2 2 5) P 1.31945, b opt

(17)

II c opt 5 As (4Ï2 2 3 1 Ï4Ï2 2 5) P 1.73366,

with maxdiff II 5 As (1 2 Ï4Ï2 2 5) P 0.09477.

P 0.10245. (20)

Note that for a ; 1 Case II is actually slightly better than Case I (maxdiff I* . maxdiff II*).

2 , z e2 , 0), E eII2 5 (2b 2 c 2 Ï1 2 2(c 2 b)2 )M for (M, y emax max

E eII4 5 (c 2 Ï3)M

P 1.80621,

(18)

II II II 1 c opt 5 2bopt , so that the solution is on the Note that a opt border of the allowed area (see Fig. 2). Note also that this

3.3. Integer Approximations Using real-valued local distances is generally not computationally desirable. Integer local distances are preferable. Candidate integer approximations of the optimal values, denoted A, B, and C, are found by multiplying the optimal local distances by an integer scale factor and rounding to the nearest integer. Then the maximal differences are computed (all expressions are available from the computations of the optimal local distances), to check the approximations. The smallest local distance, a, will act as a scale factor; therefore the resulting WDT will become k1, B/A, C/Al. The best approximation result is maxdiff* and the optimal local distances to be multiplied by the scale factor are b*opt and c*opt . Remember that maxdiff I , maxdiff II but maxdiff I* . maxdiffII*. It is thus difficult to predict in which of the two cases the best integer approximations will be found. Good integer 3 3 3 3 3 WDTs are listed in the next section. Vossepoel suggested a method for improving integer WDTs, by multiplying the resulting distance transform by a real valued scale factor (see [13]) with corrections in [1]. Verwer uses this idea as a central part of his optimizations in 3D [12]. The concept was further formalized in [5]. In the optimal WDTs the maximum negative and positive differences are equal. This is usually not the case for the integer approximations. If the WDT is rescaled so that the negative and positive errors are ‘‘balanced,’’ then the maximum error will decrease. The scale factor is computed by finding the worst positive and negative errors and com-

DIGITAL DISTANCE TRANSFORMS

TABLE 2 Integer 3 3 3 3 3 Distance Transformations Case

a

b

c

Maxdiff

PG PG PG

1 1 1

— 1 1

— — 1

1.26795 0.41421 0.73205

I II

a Iopt II a opt

bIopt II bopt

c Iopt II c opt

0.07356 0.09477

I II

1 1

bIopt * II bopt *

c Iopt * II c opt *

0.10402 0.10245

I/II I/II I II II

2 3 8 13 16

3 4 11 17 21

4 5 13 22 ~ 23 27 ~ 28

0.29289 0.11808 0.10732 0.10652 0.10296

I/II I/II I

2 3 8

3 4 11

4 5 13

0.11111 0.10008 0.10720

Scale

0.88889 0.98560 0.99991

puting the scale factor from them. For example, the k2, 3, 4l WDT, which is both Case I and Case II, has a maxdiff of 0.2929, which is not very good (the best possible value is maxdiff II* 5 0.1024). The error expressions, (10), range from E1,1 5 0 to E2,2 5 E3,2 5 0.2929. Inserting ka, Dsa, 2al and solving 2E1,1 5 E3,2 , that is, 1 2 a 5 a 2 Ï1 2 As (2a 2 a)2, gives scale factor a 5 Kl for a ‘‘balanced’’ DT and a reduced maxdiff 5 0.1111. This is a significant improvement and rather good when compared to the best possible value, maxdiff I 5 0.0736 (as a is no longer fixed at 1).

375

mal ones with a ; 1. In fact, k16, 21, 27 ~ 28l is Case II and has a better maxdiff (0.10296) than would be possible for any Case I solution (maxdiff I* 5 0.10402). As Case II WDTs have never been considered before, this is an interesting result. Finally, Table 2 lists the rescaled versions of the three first integer WDTs. Only k2, 3, 4l is significantly improved by rescaling. On the other hand, the rescaled k2, 3, 4l has almost exactly the same maxdiff as the original k3, 4, 5l, so it is doubtful if it is useful. Note that for the rescaled integer WDTs the best possible maxdiff are the ones with free a, not the ones where a ; 1. A good way to characterize a DT is the shape of its associated disc/sphere, defined as all pixels/voxels with a distance less than or equal to the radius from a single central element. In 2D, the city block and chessboard distance discs are a diamond and a square, respectively. In 3D the D6 and D26 spheres are an octahedron and a cube, respectively. The D18 sphere is a dymaxion (or cuboctahedron). In 2D, the regular 3 3 3 WDT discs are octagons. In 3D, the regular 3 3 3 3 3 WDT spheres can have any of three shapes. These are exemplified in Fig. 8. We assume that the spheres are big enough so that any jaggedness due to digitization is invisible. On top is the k3, 4, 5l sphere, which is on the border between Case I and Case II. It has 24 quadrilateral faces. At bottom left in Fig. 8 is a Case I sphere and at bottom right a Case II sphere. Here, each face of the ‘‘mixed’’ sphere has been split into two triangular faces, but the split is between different corners in the two cases. The mixed and Case I spheres can be found in previous works [12, 2], but the Case II sphere is new. All these polyhedrons are reasonably good approximations of the Euclidean sphere.

4. RESULTS

In this section the results of the optimality computations are summarized and illustrated. Table 2 lists a number of distance transforms. First, the three PGDTs are listed, with their associated maxdiff. They are easily computed from the expressions (1)–(3). The PGDTs are efficient to compute, but not good approximations of the Euclidean distance. Worse is that they are very rotation dependent. Next in Table 2 come the optimal values for Case I and Case II. The third group of results in Table 2 are the optimal values with a ; 1. After the real-valued WDTs the best integer approximations, using scale factors (5A) up to 20, are listed in Table 2. It can be noted that after scale factor 3, k3, 4, 5l, the improvement is small. Note also that in Case II, the local distance C can have several values, which reflects the fact II that c opt * can vary within an interval without affecting maxdiff (see (19)). The smallest possible maxdiffs are the opti-

FIG. 8. Three weighted distance spheres. The top one is both Case I and Case II. The bottom left is Case I and the bottom right is Case II.

376

GUNILLA BORGEFORS

5. COMPUTING DISTANCE TRANSFORMS

This section is just a quick reminder of how the PGDTs and WDTs are computed sequentially. The computations start by initializing the image. Voxels in F , from which the distances are to be computed, are set to zero and pixels in F are set to infinity, i.e., a suitably large number. The DTs are computed using only two raster scans over the initial image, independently of the size and shape of F. In the first raster scan, starting one voxel from the border, a new value is computed for each infinity-valued voxel. Compute the sums of the values of the already visited neighbors and the corresponding local distances; e.g., add b to edge neighbors. The new voxel value is the minimum of the 13 sums. This scan computes the distances from leftup-top. In the second scan, the scan direction is reversed. Again, for each voxel the sum of the neighbor’s values and the local distances are computed. The only difference is that now the voxel itself must be included, adding zero. Thus the new, final, value is the minimum of 14 values. The second scan computes the distances from right-downbottom. The algorithm in pseudo-code becomes

most previous authors that the ‘‘best’’ integer 3D weighted distance transform is k3, 4, 5l. To get even better approximations in 3D it is natural to extend the neighborhood to 5 3 5 3 5 (in 2D this led to excellent results [4]). Verwer has published some 3D results for the extended neighborhood [12]. However, there is no discussion on the regions in parameter space where his solutions are valid, which means that there is no assurance that the integer approximations he suggests are actually within the boundaries where the equations are valid. This is a problem for further study. If better accuracy than that of the k3, 4, 5l, that is, better than about 12%, is necessary, using the Euclidean distance transform may be a better choice [7, 9]. The drawback is that a vector must be stored in each pixel instead of a single value and more scans through the image are necessary for computation. Note, also, that many applications, such as skeletonizing, become much more complex using the Euclidean DT, compared to an integer-valued DT. The limited number of possible distance values usually simplifies subsequent operations considerably. REFERENCES

First pass: for i 5 2, 3, ..., lines do for j 5 2, 3, ..., columns do for k 5 2, 3, ..., layers do v 1i, j,k 5 min hv 0p,q,r 1 dnj (p,q,r)[ N

Second pass: for i 5 lines21, ..., 2, 1 do for j 5 columns21, ..., 2, 1 do for k 5 layers21, ..., 2, 1 do v 2i, j,k 5 min hv 1p,q,r 1 dnj (p,q,r)[ N

where v 0 is the initial image, N is the set of already visited neighbors in the scan, and dn is the appropriate local distance, a, b, or c. In the PGDTs the only neighbors taken into account are, of course, those from which steps can be taken in the minimal path; e.g., in D6, only the three visited area neighbors are investigated in each scan. 6. CONCLUSIONS

The geometry and equations resulting from finding minimal paths in 3D and computing their associated distance transforms when using the 3 3 3 3 3 neighborhood have been described. A new type of weighted distances has been discovered. Optimal local distances have been computed, where optimality is defined as minimizing the maximum difference from the Euclidean distance in an M 3 M 3 M image. Integer approximations have also been found. These theoretical results do not change the conclusion of

1. A. L. D. Beckers and A. W. M. Smeulders, A comment on ‘‘A Note on ‘Distance Transformations in Digital Images,’’ Comput. Vision Graphics Image Process. 47, 1989, 89–91. 2. A. L. D. Beckers and A. W. M. Smeulders, Optimization of length measurements for isotropic distance transformations in three dimensions, CVGIP: Image Understanding 55(3), 1992, 296–306. 3. G. Borgefors, Distance transformations in arbitrary dimensions, Comput. Vision Graphics Image Process. 27, 1984, 321–345. 4. G. Borgefors, Distance transformations in digital images, Comput. Vision Graphics Image Process. 34, 1986, 344–371. 5. G. Borgefors, Another comment on ‘‘A Note on ‘Distance Transformations in Digital Images,’’ CVGIP: Image Understanding 54(2), 1991, 301–306. 6. G. Borgefors, Applications of distance transforms, in Aspects of Visual Form Processing (C. Arcelli et al., Eds.), pp. 83–108, World Scientific, Singapore, 1994. 7. P.-E. Danielsson, Euclidean distance mapping, Comput. Graphics Image Process. 14, 1980, 227–248. 8. N. Okabe, J. Toriwaki, and T. Fukumura, Paths and distance functions on three-dimensional digitized pictures, Pattern Recognit. Lett. 1, 1983, 205–212. 9. I. Ragnemalm, The euclidean distance transform in arbitrary dimensions, Pattern Recognit. Lett. 14, 1993, 883–888. 10. A. Rosenfeld and J. Pfaltz, Sequential operations in digital picture processing, J. Assoc. Comput. Mach. 13(4), 1966, 471–494. 11. A. Rosenfeld and J. Pfaltz, Distance functions in digital pictures, Pattern Recognit. 1(1), 1968, 33–61. 12. B. J. H. Verwer, Local distances for distance transformations in two and three dimensions, Pattern Recognit. Lett. 12, 1991, 671–682. 13. A. M. Vossepoel, A note on ‘‘Distance Transformations in Digital Images,’’ Comput. Vision Graphics Image Process. 43, 1988, 88–97. 14. M. Yamashita and T. Ibaraki, Distance defined by neighborhood sequences, Pattern Recognit. 19, 1986, 237–246.