Precise Euclidean distance transforms in 3D from voxel coverage representation

Pattern Recognition Letters 65 (2015) 184–191

Contents lists available at ScienceDirect

Pattern Recognition Letters journal homepage: www.elsevier.com/locate/patrec

Precise Euclidean distance transforms in 3D from voxel coverage representation✩ Vladimir Ilic´ a,∗, Joakim Lindblad a, Nataša Sladoje a,b a b

Faculty of Engineering, University of Novi Sad, Trg Dositeja Obradovi´ca 6, 21000 Novi Sad, Serbia Centre for Image Analysis, Uppsala University, Box 337, 751 05 Uppsala, Sweden

a r t i c l e

i n f o

Article history: Received 28 February 2015 Available online 7 August 2015 Keywords: Distance transform Precision Coverage representation Vector propagation DT algorithm Sub-voxel accuracy

a b s t r a c t Distance transforms (DTs) are, usually, defined on a binary image as a mapping from each background element to the distance between its centre and the centre of the closest object element. However, due to discretization effects, such DTs have limited precision, including reduced rotational and translational invariance. We show in this paper that a significant improvement in performance of Euclidean DTs can be achieved if voxel coverage values are utilized and the position of an object boundary is estimated with sub-voxel precision. We propose two algorithms of linear time complexity for estimating Euclidean DT with sub-voxel precision. The evaluation confirms that both algorithms provide 4–14 times increased accuracy compared to what is achievable from a binary object representation.

1. Introduction A distance transform (DT) assigns to each image element its distance to a selected subset of image elements. Often, the selected set is an observed object in the image, and every image element is mapped to its distance to the object. The concept of DT is closely related to many procedures and operations in image processing, such as computation of morphological operations and/or geometrical representations, image segmentation, template matching, image registration, and many others. Interest for further development and improvement of methods and algorithms for DT computation, as well as for finding new applications of DTs, remains high till nowadays. Even though the main idea of DT is rather straightforward, its efficient implementation is far from trivial. Many algorithms are proposed over the years, differing in terms of time complexity, accuracy, order of pixel processing, suitability for parallelization, etc. The Euclidean distance is most often considered in DTs, due to its properties desired in many applications. The appealing property of rotational invariance of Euclidean DT (EDT) is, however, challenged by discrete binary object representations in rectangular grids, most often used in practice. We argue for efficient utilization of information contained in acquired images, as a way to address these challenges. Intensity values in original images often relate to the

✩

This paper has been recommended for acceptance by Egon L. van den Broek. Corresponding author. Tel.: +381214852280. ´ joakim@ E-mail addresses: [email protected], [email protected] (V. Ilic), cb.uu.se (J. Lindblad), [email protected] (N. Sladoje). ∗

http://dx.doi.org/10.1016/j.patrec.2015.07.035 0167-8655/© 2015 Elsevier B.V. All rights reserved.

© 2015 Elsevier B.V. All rights reserved.

(partial) coverage of each image element by the imaged object. If this information is utilized to provide coverage object representations [13], object boundary position can be estimated with sub-pixel precision. This information can be successfully utilized in EDT estimation, as was shown in [4] for the 2D case. In this paper we propose to utilize the coverage approach and volume sampled image representations to improve the performance of vector propagation based EDT algorithms in 3D. We present two algorithms for computing approximate EDT in 3D which, as shown in the tests conducted, provide between 4 and 14 times lower error of distance values compared to a binary DT, at a reasonable (linear) computational cost, (O(N) where N is the number of the image voxels). The paper is organized as follows: Section 2 provides a brief overview of related work on DTs and on the Coverage model. In Section 3 we present a novel method for estimating 3D EDT with subvoxel accuracy. Section 4 contains results of performance evaluation of the proposed method. Section 5 summarizes and concludes the paper. 2. Background 2.1. Distance transforms The distance transform (DT) computed on a grid D (commonly a subset of Zn ) w.r.t. the set S ⊂ Rn is a mapping that assigns to each grid point p ∈ D its (minimal) distance (according to a given distance metric d) to the set S:

DS ( p) = min{ d( p, q) | q ∈ S}.

(1)

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191

Traditionally, S is an object in a binary segmented image defined on D. The underlying pointwise distance d can be of different types: Manhattan, Euclidean, Chessboard, and different weighted distances. Most common are arguably DTs based on the Euclidean metric; this high interest comes primarily from its very appealing property of isotropy, leading to rotational invariance of derived shape descriptors. It is relatively straightforward to develop an exhaustive method for computing DT, but such is in most cases of too high time complexity for practical usage; depending on the image content, it is bounded by the worst case O(n4 ) and base case O(n2 ), for n being number of pixels, [3]. One way to overcome this difficulty has been to propose methods that approximate EDTs. Consequently, an important question about any proposed DT algorithm is its accuracy w.r.t. the required computational time. Among most popular algorithms are the weighted DTs proposed for Chamfer DT [1], which provide a reasonable approximation of EDT, at a very high speed. The algorithm proposed in [2] is a well-known vector propagation based approximate DT algorithm, where closest feature vectors are propagated through a sequence of raster scans. This is still a popular algorithm, providing high accuracy of distance values. Nowadays, several algorithms for computing the exact EDT of a binary image in linear time exist. A comparative survey of algorithms can be found in, e.g., [3] for the 2D case, and in [5] for the 3D case. The appealing properties of EDT (such as rotational invariance) are, however, unavoidably limited by the amount of information about shape geometry that is preserved in the object representation. In the processes of binarization and discretization a considerable amount of that information is lost; this holds in particular on the object boundary, which significantly affects DT values. One way to preserve and further exploit information present in the original images is to utilize a coverage representation [13], instead of performing image binarization. Image elements are, in this model, assigned values proportional to the portion of their volume covered by the observed object. This information can be utilized for estimating the position of the object boundary with sub-voxel precision. Work presented in [4] shows the advantages of utilizing coverage representation for improved DT computation. They present how area coverage values can be used to compute a 2D EDT which has much higher accuracy than what is achievable from a binary representation. Initial studies conducted in 3D, and presented in [7] encourage further development of DT algorithms adjusted to coverage model, and providing sub-voxel precision. Whereas [7] suggest to use discretization grids alternative to Cartesian (such as Body-Centered Cubic (BCC) and Face-Centered Cubic (FCC)) to provide higher accuracy of sub-voxel DT, we prefer to stay with the most widely used grid and to improve performance of the sub-voxel DT algorithm by increasing precision of the boundary estimation utilizing available coverage information.

Assuming a Euclidean space Rn as a reference set and a Voronoi tessellation of Rn given by a subset of Zn , where we denote the Voronoi region of a grid point x = (x1 , x2 , . . . , xn ) ∈ Zn by σ (x). Coverage digitization in Rn is defined in [13] as. Definition 1. For a given continuous object S ∈ Rn , inscribed into an integer grid Zn , the coverage digitization of S is



  (x, v(x))  x ∈ Zn ,

v(x) =

| σ (x) ∩ S| , | σ (x)|

(2)

where |A| denotes area/volume/Lebegues measure of a set A. Special cases are pixel coverage digitization for n = 2 and voxel coverage digitization for n = 3. An approximation of coverage values can be obtained by supersampling, and applying Gauss centre point digitization to the

(b)

(c)

Fig. 1. (a) Normal direction n, signed distance dv , volume coverage v (bolded volume), and distance x indicting where the object boundary intersects the x-axis. (b) Volumes v1 , v2 , v3 , v4 , v5 and v6 when ny + nz ≤ nx and (c) when ny + nz > nx .

sub-elements. Coverage value of an element is then approximated by the fraction of its sub-elements included in the Gauss centre point digitization (i.e., having centres covered by the observed object). Formal definition and more details about the coverage model can be found in [13]. Coverage digitization provides a representation of the object where spatial elements completely covered by background are assigned value 0, those completely covered by object are assigned value 1, and those appearing on the object boundary, are assigned values between 0 and 1, in accordance to their coverage by the imaged object. Several methods for estimation of coverage values in real images (i.e., coverage segmentation methods) have been proposed [6,8,11,12]. 3. Euclidean distance transforms in 3D with sub-voxel accuracy We propose novel methods for estimating 3D EDT with sub-voxel accuracy that can be used to improve any vector propagation DT algorithm. We utilize a voxel coverage representation, as described in Section 2.2. Voxel values, proportional to the relative volume of the voxel that is covered by the observed object, are used to estimate the position of an object boundary inside the voxel. 3.1. Existing estimation of sub-voxel boundary distance using voxel coverage information Under the assumption that the boundary of an object covering the voxel is locally planar (as illustrated in Fig. 1a), there exists a relation between the signed distance from the boundary to the center of the voxel dv , the volume coverage v and the normal direction n of the planar object boundary. For all normal directions n, v = 0.5 implies dv = 0 while for each normal direction parallel to one of the coordinate axes there holds a linear relation dv = 0.5 − v. For other normal directions this relation provides an approximation of the signed distance dv . This approximation was suggested in [7] as a reasonable estimation of dv when the orientation of the boundary is not known. Here, we denote this estimation by dvlin :

dvlin = 0.5 − v.

2.2. The coverage digitization model

DC (S) =

(a)

185

(3)

3.2. Novel estimation of sub-voxel boundary distance using information about voxel coverage and normal direction The position of a (planar) boundary of the object within the voxel is uniquely determined by its unit normal direction n = (nx , ny , nz ) and its signed distance dv from the voxel centre (where the sign is s.t. v < 0.5⇒dv (v) > 0). In the following we present an algorithm for computing the exact distance dv as a function of the voxel coverage v and the normal direction n of the planar object boundary. It is easy to see that, for any n, the distance dv (v) is antisymmetric in v = 0.5 and dv (0.5) = 0; it holds that dv (1 − v) = −dv (v). Hence, it is enough to observe 0 ≤ v ≤ 0.5. Also, it is sufficient to consider normal directions n = (nx , ny , nz ) such that nx ≥ ny ≥ nz ≥ 0. The other cases are based on symmetry and can be derived by changing the sign of nx and/or ny and/or nz or if necessary by their swapping.

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191

By translating a plane with normal direction n = (nx , ny , nz ), from the position in which it intersects the vertex at (0, 0, 0) to the position where it intersects the vertex at (1, 1, 1), we determine, in order, the volumes v1 , v2 , v3 , v4 , v5 and v6 (shown in Fig. 1) as the parts of the voxel to the left of the plane when it intersects the respective voxel vertex. The order in which the plane passes through the different vertices of the voxel depends on the relation between nx and ny + nz . In the case when ny + nz ≤ nx this order is (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1) and (1, 1, 0) (Fig. 1b), while in the case when ny + nz > nx , the order becomes (0, 0, 1), (0, 1, 0), (1, 0, 0), (0, 1, 1), (1, 0, 1) and (1, 1, 0) (Fig. 1c). For a given normal direction n it holds that:

v1 =

1 n2z , 6 nx ny

v2 =

1 n2z − 3nx nz + 3n2x 1 (nx − ny )3 − , 6 nx ny 6 nx ny nz

v4 = 1 − v3 ,

v5 = 1 − v2 ,

v6 = 1 − v1 .

(4)

Now we can present our main contribution of the paper: a complete algorithm which, for a given normal direction n, computes the signed distance dv , as a function of a volume coverage v, 0 ≤ v ≤ 0.5. Algorithm 1.

if (ny = 0) /* implies nz = 0 */ dv = 0.5 − v elseif (v ≤ v1 )  dv = 12 (nx + ny + nz ) − 3 6v nx ny nz elseif (v ≤ v2 )  1 2 (nx

n3z n3y

c=3

+ ny ) −

b = −3

;

n3z (n2y +n2z ) n2x n3y

2 v nx ny −

1 2 12 nz

;

(ny +nz )n3z nx n3y

d = 6v

n4z n2x n2y

x = CardanRoots(a, b, c, d) if (ny + nz ≤ nx) n n +n select x ∈ [ nyx , ynx z ] else n select x ∈ [ nyx , 1] endif dv = 12 (nx + ny + nz ) − xnx else if (ny + nz ≤ nx) n +n x = v + 12 ynx z else n3 n3 b = −3 3z a = 2 z3 ; ny

c=3

0 −0.2 −0.4 −0.6 −0.8 −1 0

0.1

2 2 2 n3z nx +ny +nz n2x n3y

ny

;

−

n3z (n3y +n3z ) n3x n3y

0.2

0.3

0.4 0.5 0.6 Volume coverage v

0.7

0.8

0.9

1

b , 3a √ b i 3 S+T x2 = − − + (S − T ), 2 3a 2 √ b i 3 S+T − − (S − T ), x3 = − 2 3a 2 x1 = S + T −



S=

3

R+



 Q3

+

R2 ,

T=

3

R−



Q 3 + R2 ,

and

Q=

elseif (v ≤ v3 and nz = 0) a=

0.2

where

Input: normal direction of the boundary, n = (nx , ny , nz ) with nx ≥ ny ≥ nz , and voxel coverage v, where 0 ≤ v ≤ 0.5. Output: signed distance dv from the voxel centre to the object boundary.

dv =

0.4

Computation of dv by Algorithm 1 requires, at two instances, to compute the roots of a third degree polynomial. We do this by calling a function CardanRoots(a, b, c, d), which solves a cubic equation ax3 + bx2 + cx + d = 0 for a, b, c, d ∈ R utilizing Cardano’s formula. The roots, x1 , x2 , and x3 of this equation are given by

1 ny + nz , 2 nx

else

v3 =

n=(1,0,0)/|(1,0,0)| n=(1,1,0)/|(1,1,0)| n=(1,1,1)/|(1,1,1)|

0.6

Fig. 2. Signed distance function dv (v), 0 ≤ v ≤ 1, for three different unit normal directions n = (nx , ny , nz ).

1 n2z − 3ny nz + 3n2y , 6 nx ny

if (ny + nz ≤ nx )

v3 =

1 0.8

Signed distance dv

186

3ac − b2 , 9a2

R=

9abc − 27a2 d − 2b3 . 54a3

How to choose the root appropriate in our task is, for each of the relevant cases, explicitly stated in Algorithm 1. Plots of the function dv (v), 0 ≤ v ≤ 1, for three given unit normal directions n = (nx , ny , nz ), are shown in Fig. 2. As observed, function dv (v), 0 ≤ v ≤ 1 is continuous, monotonically decreasing and such that dv (1 − v) = −dv (v). We can also notice that dv can be fairly well approximated by a linear function dvlin = 0.5 − v which is equivalent to dv for normal direction n = (1, 0, 0). A more complete visualization of the function dv (n, v) is shown in Fig. 3 where distances dv for 6 different voxel coverage values are plotted as functions of √ normal direction. The top surface, corresponding to v = 0 peaks at 3/2 and the bottom one, for v = 0.5 is constant equal to zero. The complete function is smooth and can be well tabulated and interpolated for increased speed. 3.3. Proposed 3D Euclidean distance transforms

nx +ny +nz nx

d = 6v

x = CardanRoots(a, b, c, d) n +n select x ∈ [1, ynx z ] endif dv = 12 (nx + ny + nz ) − xnx endif

n4z n2x n2y

−

3 3 3 n3z nx +ny +nz n3y n3x

The signed distance dv from the center of the voxel to the boundary of the observed object within that voxel, as computed by Algorithm 1, can now be used to improve the Euclidean distance estimation of any vector propagation based EDT. We can, therefore, formulate a novel method for EDT estimation, with sub-voxel accuracy. The classic binary vector propagation EDT computes for each background voxel a vector of integers d i = (xi , yi , zi ) pointing to the closest object voxel. A binary EDT (BEDT) assignes to that voxel, the value computed as:

di =



x2i + y2i + zi2 .

(5)

An approach to improve the accuracy of BEDT is suggested in [7]. The linear approximation, dvlin , given in (3), of the distance dv is used to

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191

187

0.9 0.8 0.7

distance dv

0.6 0.5 0.4 0.3 0.2 0.1 pi/2

0 0

pi/4

pi/4

pi/2

0

polar angle θ

azimuth φ

Fig. 3. Signed distance function dv (n, v) for v ∈ {0, 0.1, 0.2, 0.3, 0.4, 0.5} (top to bottom) and normal directions in the first octant, parametrized by spherical coordinate angles  and .

define AAEDT (Anti-aliased Euclidean Distance Transform):

dAAEDT = dVCEDT 1 = dBEDT + dvlin = dBEDT + 0.5 − v.

(6)

For consistency, we will, in this paper, refer to this approximation as VCEDT1 (Voxel Coverage EDT 1). This approach utilizes only the local information about the intensity of the particular voxel, and is therefore simple and reasonably fast. It is clear that information about the normal direction in the boundary voxels can further improve the accuracy of 3D EDT. That is the approach we take in this paper. One possibility to incorporate the information of a normal direction into the distance computation is to assume that the boundary of the object is approximately orthogonal to the distance vector di assigned to that voxel. This assumption is asymptotically correct; it may be violated for voxels close to the object boundary, whereas it provides a good approximation of the boundary position for voxels further away. Utilizing this assumption, we compute dv by Algorithm 1, and denote it dvort . We then, in computed EDT, assign distance value d to each background voxel as

dVCEDT 2 = dBEDT + dvort .

(7)

Further improvement can be made by taking additional care of the distance computation close to the object boundary. We estimate the unit normal direction n in each boundary voxel by using three 3 × 3 × 3 Zucker–Hummel gradient filters [14]. The sub-voxel distance dv computed by Algorithm 1 using such a gradient based normal direcgrad tion estimate is denoted dv . For voxels close to the object boundary grad we use dv to improve the EDT estimate. For voxels further from the boundary we rely on the same approach as in dVCEDT2 . Denoting

db = |di × n | ,

da = di · n ,

dVCEDT 3 =

da + dvgrad , dBEDT + dvort ,

for |dvgrad | ≤ 0.5 ∧ db ≤ 0.5 otherwise.

with a separable (squared) EDT, such as [10] or [9] is not equally straightforward. In the following we evaluate performance of our proposed methods dVCEDT2 and dVCEDT3 , in comparison with the classic dBEDT , as well as dVCEDT1 [7]. 4. Evaluation Our test set consists of Gauss centre point respective voxel coverage digitizations of spheres with 30 different real-valued diameters ranging from 2 to 62 voxels and cubes with real-valued edge lengths ranging from 2 to 42 voxels. For each diameter, and edge length, we generate, respectively, 50 spheres, and 50 cubes, with centres randomly positioned within the voxel; cubes are observed in random orientation (a composition of 3 successive random rotations about the coordinate axes). The objects are digitized in a rectangular grid of a size 80 × 80 × 80 voxels. Volume coverage of a voxel is estimated utilizing super-sampling of boundary voxels by a factor 16, and counting the number of sub-voxel centers covered by the object (as explained in details in [13]). We compare the performance of our suggested methods both with the classical sequential EDT vector propagation algorithm using 26neighbors, presented in [1], and with the approach of [7], where computation of dVCEDT1 utilizes information about coverage but not of normal direction. We compute DTs of our test objects by the four methods (referred to as BEDT, VCEDT1, VCEDT2, VCEDT3). The same propagation code is used for BEDT and voxel coverage methods, with the additional improvement of the later by the considered sub-voxel distance, as described above. Three quantitative performance measures are observed for each of the methods: the root-mean-square error (RMSE)



RMSE =

(8)

n 1 (xˆi − xi )2 , n

(10)

i=1

the mean-absolute-error (MAE)

(see Fig. 4) we define dVCEDT3 as



Fig. 4. Local gradient vector g, signed sub-voxel distance dvgrad , distance di (Eq. (5)) and distances da , and db (Eq. (8)).

(9)

Both proposed algorithms are of linear complexity w.r.t. the number of voxels in the image. Algorithm 1, which is evaluated for each voxel, does not contain any loops and has constant O(1) complexity. To estimate the object normal direction we utilize the distance vector di , which must be available when computing the distance value. In the vector propagation algorithm the distance vector is directly available at each step. To combine our proposed framework

MAE =

n 1 |xˆi − xi | , n

(11)

i=1

and the empirical range of errors

Range = max (xˆi − xi ) − min (xˆi − xi ) , i=1...n

i=1...n

(12)

where n is size of the image, xi are the true values of EDT and xˆi are the estimated EDT values. The computed RMSE, MAE and Range values are presented in Fig. 5. For each size (diameter, or edge length) 50 spheres and 50

188

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191

0.7

0.7 BEDT VCEDT1 VCEDT2 VCEDT3

0.6 0.5

0.5

0.4 0.3

0.4 0.3

0.2

0.2

0.1

0.1

0 0

10

20

30 40 Diameter

50

60

BEDT VCEDT1 VCEDT2 VCEDT3

0.6

RMSE

0 0

70

10

(a) RMSE

40

0.5

0.5 MAE

0.3

RMSE (voxels) MAE (voxels) Range (voxels) Rel. improvement Rel. improvement

0.2016 0.1371 0.9999 1 –

0.1054 0.0718 0.5994 1.91 1

0.0216 0.0115 0.5063 11.92 6.24

0.0181 0.0101 0.3042 13.57 7.11

Table 2 RMSE, MAE and Range (measured in voxels), and relative improvement (w.r.t. MAE) of BEDT, VCEDT1 and proposed VCEDT2, VCEDT3, applied to a cube with edge length 51.0938.

50

0.3

0.2

0.2

0.1

0.1 10

20

30 40 Diameter

50

60

0 0

70

10

20 30 Edge length

40

50

(d) MAE

1.4

1.4

BEDT VCEDT1

1.2

1.2

VCEDT2 VCEDT3

1

1

0.8

0.8

Range

Range

VCEDT3

0.4

(c) MAE

0.6

0.6

0.4

0.4

0.2

0.2

0 0

10

20

30 40 Diameter

50

60

BEDT VCEDT1 VCEDT2 VCEDT3

0 0

70

10

(e) Range

20 30 Edge length

40

50

(f) Range

Fig. 5. RMSE, MAE and Range of errors, measured in voxels, of BEDT, VCEDT1 and the proposed VCEDT2, VCEDT3 applied to spheres (first column) and cubes (second column) of increasing sizes.

5

4

4

VCEDT2

VCEDT3

RMSE (voxels) MAE (voxels) Range (voxels) Rel. improvement Rel. improvement

0.2851 0.2048 1.3140 1 –

0.1500 0.1173 0.9958 1.75 1

0.0954 0.0515 0.9813 3.98 2.28

0.0946 0.0503 0.7774 4.07 2.33

5

x 10

4

2.5 2 1.5

3

2.5 2 1.5

2.5 2 1.5

1

1

1

1

0.5

0.5

0.5

0.5

0.2

0.4 0.6 Errors

0.8

0 −0.4

1

−0.2

(a) BEDT

0.2

0 −0.4

0.4

3

x 10

3

2

1

Frequency

2.5

2

Frequency

2.5

2

0.5

1.5 1 0.5

0.5

1 Errors

(e) BEDT

1.5

0 −0.4

0.2

0

0.2 Errors

0.4

(f) VCEDT1

−0.2

0.6

0.8

0.2

0.4

(d) VCEDT3

x 10

3

x 10

2.5

1.5 1

0 −0.4

0 Errors

5

0.5

−0.2

0 −0.4

0.4

5

2.5

1.5

0 Errors

(c) VCEDT2

5

x 10

0 0

−0.2

(b) VCEDT1

5

3

0 Errors

Frequency

0 0

x 10

3.5

3

Frequency

2 1.5

VCEDT1

3.5

3

Frequency

3

BEDT

5

x 10

3.5

2.5

Methods

cubes, each with a random orientation and placement in the digital grid, are observed. The errors, measured in voxels, are averaged over these 50 instances. It is clear that the use of volume coverage values leads to improvement in accuracy over the BEDT. Our proposed VCEDT2 and VCEDT3 in general outperform VCEDT1, particularly in the case of spheres. This improvement is present for larger cubes as well, whereas the proposed methods do not provide improvement for small cubes (with edge length less than 20 voxels). The reduction in range of errors achieved by method VCDTE3 is clearly present for all the observed objects. For further empirical evaluation we observe the errors of the computed EDT of a randomly positioned sphere with a real-valued diameter of 81.8268 voxels, and a randomly positioned and rotated cube with real-valued edge length of 51.0938 voxels. Rectangular sampling grid of a size 100 × 100 × 100 voxels is considered. We show histograms of errors for these objects in Fig. 6 and provide detailed results in Table 1 (for the sphere), and Table 2 (for the cube). Significant increase of both precision and accuracy, in case of VCEDT2 and VCEDT3, compared to both BEDT and VCEDT1, is clearly visible. The

5

x 10

3.5

Frequency

VCEDT2

BEDT VCEDT1 VCEDT2 VCEDT3

0.6

0.4

0 0

Frequency

VCEDT1

0.7 BEDT VCEDT1 VCEDT2 VCEDT3

0.6

4

BEDT

(b) RMSE

0.7

MAE

20 30 Edge length

Methods

Frequency

RMSE

Table 1 RMSE, MAE and Range (measured in voxels), and relative improvement (w.r.t. MAE) of BEDT, VCEDT1 and proposed VCEDT2, VCEDT3, applied to a sphere of diameter 81.8268.

2 1.5 1 0.5

−0.2

0

0.2 Errors

0.4

(g) VCEDT2

0.6

0.8

0 −0.4

−0.2

0

0.2 Errors

0.4

0.6

0.8

(h) VCEDT3

Fig. 6. Histograms of errors of BEDT, VCEDT1 and proposed VCEDT2, VCEDT3, when applied to a sphere of diameter 81.8268 (first row) and a cube with edge length 51.0938 (second row).

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191

189

Fig. 7. First row: the 2D cross-sections of 3D error-image, with the largest sum of absolute errors, for (a) BEDT, (b) VCEDT1, (c) VCEDT2, (d) VCEDT3, applied to a sphere (diameter 81.8268). Second row: histograms for respective 2D cross-sections shown in the first row.

0.35 BEDT VCEDT1 VCEDT2 VCEDT3

0.3

0.25

MAE

0.2

0.15

0.1

0.05

0 0

20

40

60

80

100

120

140

160

180

Angle

Fig. 9. Al Capone object. Fig. 8. MAE of BEDT, VCEDT1 and proposed VCEDT2, VCEDT3, when applied to a randomly positioned cube with edge length 50.6802 which is rotated around the z-axis from 0 to 180 degrees with a step of 5 degrees.

Table 3 RMSE, Range, MAE (in voxels), of the cross-section of errorimage with the largest sum of absolute errors, for BEDT, VCEDT1 and proposed VCEDT2, VCEDT3, applied to a sphere. Methods

BEDT

VCEDT1

VCEDT2

VCEDT3

RMSE MAE Range

0.2282 0.1890 0.8180

0.1335 0.1089 0.4482

0.0210 0.0139 0.2254

0.0194 0.0131 0.2254

results indicate 4–14 times increased accuracy (reduction in MAE) by utilizing the proposed methods VCEDT2, VCEDT3, compared to BEDT, and 2–7 times increased accuracy relative to VCEDT1. For the sphere, 46.72%, 66.94%, 99.28% and 99.62% of the voxels are assigned distance values which are accurate up to ±0.1 voxels, in the BEDT, VCEDT1, VCEDT2 and VCEDT3 distance maps, respectively. For the cube, we observe that 61.98%, 78.05%, 93.11% and 93.25% of the voxels are assigned distance values with deviations up to ±0.2 voxels, in the BEDT, VCEDT1, VCEDT2 and VCEDT3 distance maps, respectively. Performance of the observed methods is further visualized in Fig. 7, where the errors of a 2D slice of the generated 3D distance

Table 4 MAE (in voxels) of BEDT, VCEDT1, VCEDT2, and VCEDT3, when applied to a randomly positioned cube with edge length of 50.6802 voxels which is rotated (top rows) around the z-axis from 0 to 180 degrees with a step of 5 degrees or translated (bottom rows) to 20 random positions. Methods

BEDT

VCEDT1

VCEDT2

VCEDT3

Rot. MAE (average) Rot. MAE (std. dev.) Transl. MAE (average) Transl. MAE (std. dev.)

0.2003 0.0142 0.1955 0.0066

0.1065 0.0095 0.0941 0.0014

0.0533 0.0022 0.0574 0.0015

0.0522 0.0023 0.0563 0.0015

maps for one sphere are shown. We show in the first row of Fig. 7, for each method, the 2D section for which the sum of absolute errors of the distance values assigned to the voxels is the highest. For all distances, the worst case error 2D section corresponds to a plane close to the object boundary. Corresponding histograms of deviations are presented in the second row of Fig. 7. Table 3 presents values of the observed errors computed for this particular 2D section. We notice a clear decrease of errors when the proposed methods are used to estimate EDT. To further emphasize the improved rotational and translational invariance achieved with the proposed methods, we observe the MAE of BEDT, VCEDT1 and proposed VCEDT2, VCEDT3, when applied to a randomly positioned and rotated cube with a real-valued edge length of 50.6802 voxels in a sampling grid of a size 100 × 100 × 100 voxels.

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191 14000

14000

12000

12000

12000

10000

10000

10000

10000

8000

6000

8000

6000

8000

6000

4000

4000

4000

2000

2000

2000

0 -0.5

0

0.5

1

1.5

2

0 -0.4

2.5

-0.2

0

Errors

0.2

0.4

0.6

0.8

Frequency

14000

12000

Frequency

14000

Frequency

Frequency

190

0 -0.4

6000

4000

2000

-0.2

0

Errors

(a) BEDT

8000

0.2

0.4

0.6

0.8

0 -0.4

-0.2

0

Errors

(b) VCEDT1

(c) VCEDT2

0.2

0.4

0.6

0.8

Errors

(d) VCEDT3

Fig. 10. Histograms of errors for BEDT, VCEDT1, VCEDT2, and VCEDT3 when computed for 50 × 50 × 50 Al Capone object.

shown execution time vs. image size. Further increase of speed can be achieved by pre-computing a lookup table for the distance function computed by Algorithm 1.

45

40

Execution time (seconds)

35

30

5. Conclusion BEDT VCEDT1 VCEDT2 VCEDT3

25

20

15

10

5

0 0

0.5

1

1.5

2

2.5

3

3.5 × 106

Number of voxels

Fig. 11. Execution time (in seconds) of BEDT, VCEDT1, VCDET2 and VCEDT3, when applied to Al Capone object of size 183 , 303 , 503 , 903 , and 1503 voxels, plotted as a function of the number of voxels (N) in the image. Table 5 RMSE, MAE and Range (measured in voxels), of BEDT, VCEDT1, VCEDT2, and VCEDT3, when computed for 50 × 50 × 50 Al Capone object. Methods

BEDT

VCEDT1

VCEDT2

VCEDT3

RMSE MAE Range

0.1858 0.3510 2.1773

0.0124 0.0819 1.0068

0.0119 0.0771 1.0128

0.0118 0.0760 0.9037

This cube is then additionally rotated around the z-axis for angles from 0 to 180 degrees with a step of 5 degrees. We show MAE for, in this way, generated cubes as a function of the rotation angle in Fig. 8. The absolute reduction of the error, as well as the reduced angular variation, is clearly visible. We also performed test with cubes randomly displaced within one voxel, with similar results, with the exception that VCEDT1 has equally low translational variance as VCEDT2 and VCEDT3 (still at a higher overall error-level). Table 4 presents average MAE over all rotations and translations, as well as standard deviation of MAE, showing how the distance measure varies with rotation and translation, respectively. A significant decrease of both rotational and translational variability, of proposed VCEDT2 and VCEDT3, compared to both BEDT and VCEDT1, is clearly visible. To evaluate the performance on a more complex shape, we created, by 9-times sub-sampling, a 50 × 50 × 50 voxel coverage representation of the Al Capone object (Fig. 9) from the IAPR-TC18 webpage (http://www.tc18.org), courtesy of E. Remy. The original 450 × 450 × 450 voxel object is used to compute the “ground truth” distance values using the method of [9]. Distribution of distance errors for BEDT, VCEDT1, VCEDT2, and VCEDT3 are shown in Fig. 10, and error measures are listed in Table 5. Increase in accuracy from use of coverage information is clear, while the improvement over VCEDT1 is less prominent than for simpler objects. The proposed algorithm is of linear time complexity w.r.t. number of voxels in the image. We illustrate this property in Fig. 11, where is

We present two novel 3D EDT estimations methods of linear time complexity which utilize voxel coverage information to achieve increased accuracy and precision of the estimated distance values. Voxel coverage values, which are proportional to the relative volume of the voxel covered by the observed continuous object, are used for estimating the position of the boundary of the object inside the voxels, with sub-voxel precision. This enables significantly improved accuracy of estimated distance values assigned in EDT. We have performed a statistical evaluation of the proposed methods, comparing them with both classical binary vector propagation based EDT, and with an existing method that relies on voxel coverage representation and considers sub-voxel position of the object boundary. The conducted tests show clearly the improvement in accuracy and precision when our proposed methods are applied. The proposed VCEDT3 slightly outperforms VCEDT2 by reducing the worst errors and thus shrinking the range of errors. Both VCEDT2 and VCEDT3 outperform both BEDT and VCEDT1; mean absolute error is reduced by a factor up to 14 compared to BEDT, and up to 7 compared to VCEDT1. We conclude that the advantage of utilizing coverage information in images is confirmed one more time. An appropriate treatment of that information, through the suggested EDT estimation methods, leads to increased accuracy and precision of the distance values. Our further work will be devoted to adjustment and development of image analysis methods that can benefit from this well performing EDT estimator as well as further evaluation of the proposed methods, when applied on images resulting from real imaging conditions. Acknowledgments Authors are supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia. V. Ilic´ is supported through Project ON174019. J. Lindblad and N. Sladoje are supported through Projects ON174008 and III44006. N. Sladoje is supported by the Swedish Governmental Agency for Innovation Systems (VINNOVA Grant No. 2014-01432). References [1] G. Borgefors, Distance transformations in arbitrary dimensions, Comput. Vis. Graph. Image Process. 27 (1984) 321–345. [2] P.E. Danielsson, Euclidean distance mapping, Comput. Graph. Image Process. 14 (1980) 227–248. [3] R. Fabbri, L.D.F. Costa, J.C. Torelli, O.M. Bruno, 2D Euclidean distance transform algorithms: a comparative survey, ACM Comput. Surv. 40 (1) (2008) 2. [4] S. Gustavson, R. Strand, Anti-aliased Euclidean distance transform, Pattern Recognit. Lett. 32 (2011) 252–257. [5] M.W. Jones, J.A. Baerentzen, M. Sramek, 3D distance fields: a survey of techniques and applications, IEEE Trans. Vis. Comput. Graph. 12 (4) (2006) 581–599. [6] J. Lindblad, N. Sladoje, Coverage segmentation based on linear unmixing and minimization of perimeter and boundary thickness, Pattern Recognit. Lett. 33 (6) (2012) 728–738.

V. Ili´c et al. / Pattern Recognition Letters 65 (2015) 184–191 [7] E. Linnér, R. Strand, Anti-aliased Euclidean distance transform on 3D sampling lattices, in: Proceedings of the 18th International Conference on Discrete Geometry for Computer Imagery – DGCI 2014, Siena, Italy, in: Volume 8668 of Lecture Notes in Computer Science, Springer International Publishing, Switzerland, 2014, pp. 88–98. [8] F. Malmberg, J. Lindblad, N. Sladoje, I. Nyström, A graph-based framework for subpixel image segmentation, Theor. Comput. Sci. 412 (15) (2011) 1338–1349. [9] C. Maurer, R. Qi, V. Raghavan, A linear time algorithm for computing exact Euclidean distance transform of binary images in arbitrary dimensions, IEEE Trans. Pattern Anal. Mach. Intell. 25 (2) (2003) 265–270. [10] A. Meijster, J.B. Roerdink, W.H. Hesselink, A general algorithm for computing distance transforms in linear time, in: Proceedings of Mathematical Morphology and its Applications to Signal and Image Processing, Computer, Imaging and Vision, vol. 18, Springer, 2000, pp. 331–340.

191

[11] N. Sladoje, J. Lindblad, High-precision boundary length estimation by utilizing gray-level information, IEEE Trans. Pattern Anal. Mach. Intell. 31 (2) (2009) 357– 363. [12] N. Sladoje, J. Lindblad, Pixel coverage segmentation for improved feature estimation, in: Proceedings of International Conference on Image Analysis and Processing, ICIAP, in: Volume 5716 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Heidelberg, 2009, pp. 929–938. [13] N. Sladoje, J. Lindblad, The coverage model and its use in image processing, in: Selected Topics on Image Processing and Cryptology, Zbornik radova (Collection of Papers), vol. 15, Mathematical Institute of the Serbian Academy of Sciences and Arts, Belgrade, 2012, pp. 39–117. [14] S.W. Zucker, R.A. Hummel, A three-dimensional edge operator, IEEE Trans. Pattern Anal. Mach. Intell. 3 (1981) 324–331.