Euclidian Distance Transform Algorithms
A Comparative Study
Ola Nilsson∗ Andreas S¨oderstr¨om †
Graphics Group
Link¨oping University
Abstract
We compare the most frequently used algorithms for computing distance trans-forms in terms of speed, scalability and accuracy. The algorithms we consider are: partial differential equation based distancing methods of different finite difference ac-curacy, fast marching methods and fast sweeping methods. Our results show clearly that all the first order algorithms perform similarly in terms of accuracy and that the only major difference is efficiency. In parallel with the revision of the preprint we find that the newly proposed O(N) quantized fast marching algorithm in practice does not scale linearly.
1
Theory
Implicit methods has gained wide acceptance in many applications in image and geom-etry processing. Segmentation in two and three dimensions are good examples, which essentially performs tracking of a deforming interface. The capturing of such an inter-face is often done with a so-called level-set which maps its zero level set around the interface and then moves with it. The full theory of level-sets is beyond the scope of this short paper, and we point the reader to the excellent introductions [7, 10] and the original publication [8]. For some examples of level-set applications we can recommend [6, 11]. What we will discuss here are some of the numerical issues related to the under-lying scalar function that implicitly samples the level-set function φ(~x), and a number of algorithms to keep it numerically stable.
The usual choice of scalar function is a Euclidian signed distance transform. Such a
∗
olani@itn.liu.se
†
function describes the Euclidian distance (d) as |d(~x)| = min(|~x − ~xi|)
|~a| =qa2
x+ a2y + a2z and is a signed distance in the sense that
d(~x) < 0, ~x ∈ Ω− d(~x) > 0, ~x ∈ Ω+ d(~x) = 0, ~x ∈ Γ
The domain is divided into outside, Ω+, inside, Ω−, and interface, Γ. This means that the geometry needs to form a closed surface, or be orientable, for an example see figure 1.
Ω Γ
Ω+
−
Figure 1: An example domain in <2 with the shape defined by Γ (left), and a color coded signed distance function (right).
The signed distance function is most often defined implicitly by its gradient
|∇φ(~x)| = 1, ~x ∈ <n (1)
which classifies it as a specialized Eikonal equation. It has the boundary condition φ(~x) = 0, ~x ∈ Γ ⊂ <n
In this paper we study different methods that given a boundary condition (i.e. an interface) construct a signed Euclidian distance transform by solving equation (1) over a rectangular domain in <2.
2
Method
2.1 The fast marching method
Sethians distancing method [9] is based on the upwinding nature of equation 1. This feature has the consequence that the solution at any point is only depending on points
closer to the interface. Using a data structure that exploit this it should be possible to solve the equation in linear time1 visiting each grid cell a constant number of times as the solution is only depending on neighbors already solved for.
The fast marching method (FMM) uses a heap to achieve this. First loop over the grid cell neighbors to the interface, update their values and insert them into the heap. Then solve for the smallest unsolved grid cell by
• extracting the top (smallest) element in the heap and freeze the value • update and (possibly) insert its neighbors.
This is continued until the heap is empty. The given solution is valid since everytime an element is extracted it is guaranteed that its solution is only depending on already solved for grid cells which are closer to the interface. Although the idea is linear in the number of grid cells, N , the heap has an O(log N ) insertion and update overhead making the algorithm’s total runtime O(N log N )
2.2 The quantized fast marching method
A complimentary fast marching implementation has recently been proposed by Yatziv et al. [4]. It’s based on the FMM but alleviates the need for an expensive heap by using what they call an untidy heap, which quantizes the sorting and thus goes from the insert/update in O(log N ) to O(1) making the algorithm theoretically O(N ), where N is the number of grid cells. The untidy heap is in fact a circular array of pointers to bins, implemented as linked lists, as was proposed in [2]. The circular nature of the data structure drastically reduces memory usage, and is justifiable considering the scope of values residing at any one time within the heap.
The untidy heap of course introduces errors in the algorithm, but the quantization is only used when placing elements into the bins, not for the actual computations. This limits the error of each drawn top element to be within the range of its bin. And if the number of bins is chosen in an intelligent manner the error can be bound within the same scaling as the rest of the algorithm.
In our implementation the number of bins, b, denote the number of concurrent bins meaning that at any time, the quantized heap in 2D contains only values in the inter-val
[th, th + √
2h]
divided in b bins, where t denotes how far the solution has reached and h is the grid cell size.
1
2.3 The fast sweeping method
Another way of solving equation (1) is to take advantage of the fact that although the solution at a grid cell is only depending on neighbors with smaller values the number of neighbors is limited. With four-connectivity the solution at a grid cell can only be depending on four different possible combinations of the neighbors (figure 2). This together with the fact that the solution must be the closest distance to the interface makes it possible to solve for all four stencil directions and then take the minimum of the potential solutions. Since the solution also must be able to propagate away from the
Figure 2: The four possible fast sweeping stencils.
interface it is necessary to run the method over the grid in an consistent manner. One such run has a corresponding fixed stencil and is called a sweep. To solve for a complete domain simply sweep it over and solve for each grid cell using the corresponding stencil. In 2D this entails four sweeps (see figure 3), in 3D eight, and in N D N2 sweeps. This
Figure 3: The four sweeping directions. x+y+, x-y+, x-y-, x+y-.
fast sweeping method (FSM) bears remarkable similarities to early work at Link¨oping University by Danielsson [3] and was introduced as the fast sweeping method (FSM) to the level set community by [12].
2.4 The PDE based distancing method
The idea behind this method is to rewrite equation (1) as a time-dependent PDE on the form
dφ(~x)
dτ + sign (φ(~x)) · (|∇φ(~x)| − 1) = 0 (2)
with the interface as boundary condition. As the fictious time τ increases this equation will converge towards a stable state where |∇φ(~x)| = 1. Implementing a numerical scheme that does this is straight-forward, but there are a few pitfalls. Numerically equation (2) will propagate information outwards from the boundary, i.e. our interface. To achieve this we will need to discretize our spatial differentials in such a way that information in the direction of the interface is prioritized. To do this we use Godunov’s upwinding scheme [7]. Also, to improve the behaviour of the algorithm a smeared-out
sign function can be used. We use
sign(φ(~x)) = φ(~x)
pφ(~x) + |∇φ(~x)| (3)
In this paper we will use two versions of PDE based distancing, one that is first order accurate and one that is at least 3rd order and ideally 5th order accurate.
2.4.1 First order accurate PDE distancing
For the first order accurate method we use euler integration
φ(~x)τ +∆τ = φ(~x)τ+ ∆τ f (4)
f = sign(φ(~x)) · (1 − |∇φ(~x)|) (5)
to propagate the distance field one timestep ∆τ forward in time. In order for the time integration scheme to be stable the CFL condition, ∆τ < 1, must be satisfied[7]. We use ∆τ = 0.7 which means that we will propagate information at a speed of 0.7 grid cells per iteration of the algorithm. For the spatial differentials we use
Di+= df dx ≈ fi+1− fi ∆x Di−= df dx ≈ fi− fi−1 ∆x which are first order accurate directional differentials.
2.4.2 High order accurate PDE distancing
The main advantage of the PDE based distancing algorithm is the possibility to achive highly accurate distance fields. To do this we will need to use high order accurate schemes for discretizing the spatial differentials. The scheme we use is the Hamilton-Jacobi Weighted Essentially Non-Oscillating scheme [5], or HJ WENO for short. This scheme uses a stencil of six grid-points and weighs the contribution from each in an intelligent manner to achieve at least 3rd order accuracy. Under ideal conditions the scheme will be 5th order accurate. The details of this scheme is rather complex, but once implemented it can be used together with the same algorithm described above. Our experience was however that this scheme is rather sensetive to the size of the time-step. Good results were achieved using ∆τ = 0.1.
3
Experiments
When comparing the different distancing algorithms it is important to eliminate any external sources of errors. To achieve this we decided to use a simple setup. We use
a rectangular two dimensional domain with the physical dimensions 100x100 distance units. In the center of this domain we inscribe the level set of a circle with radius 25 distance units. We then run the different distancing methods over the domain and obtain the results.
Numerically the domain is represented as a grid which is N × N . Only grid cells for which
φ(i, j) ∈ Γ explicitly indicated by
|φ(i, j)| ≤ 0.5 d.u.
(d.u. = distance units) are known from initialization, the rest of the grid is defined as plus/minus a threshold value denoting inside/outside.
3.1 Initial condition
As initial condition for the re-distancing algorithms we use the grid cells that contains the interface of the circle. These pixels store the analytically correct distance, thus we have a perfect initial condition. When doing re-initialization of distance fields of level sets this is generally not the case and one will need to initialize this closest band of pixels. This can be troublesome if high order methods are used and is not covered in the benchmark study.
3.2 Benchmarking
We will benchmark each method by computing distances for a full image initialized as described above. By successively increasing the resolution from 100x100 grid points to 1000x1000 grid points we will compare accuracy, speed and scalability of the different methods. Accuracy is measured by comparing the distanced image to the analythical distance from the circle and also by measuring the mean square error of the result. Speed is measured as run times for the different methods on the same hardware setup, and scalability is measured as time per pixel with respect to increasing grid sizes. Since the PDE based algorithms do not need to lock down the distance values of the initial grid cells we will test these methods using both a stationary and non-stationary interface.
4
Results
The results from the accuracy tests can be seen in table 1. The table contains the mean square error of the different algorithms. The distribution of the error is visualized by plotting the difference
between the correct image and the calculated image for each method. These plots can be seen in figures 5 trough 13 where the z-axis represents the error. The color coding is equalized over [min − max] for each plot. As the quantized FMM potentially is as good as the FMM special comparisons were performed to verify this. The difference between quantized FMM (QFMM) and the normal FMM is provided in figure 14, and the accuracy with respect to the number of bins in the QFMM is shown in figure 13. Run times are provided for all the methods in table 2. All timings were performed on a Macintosh G5 2.5 GHz machine with 2GB RAM2. Comparisons with respect to amortized per-pixel run times can be seen in figure 4 and show clearly how the different methods scale with respect to increased grid size. We included the quantized FMM method with both 4, and 8 bins to show how the quantization levels impacts on run time. For a discussion of the results and how these methods perform we refer to section 5.
Table 1: Mean square error in distance units
Method 100 × 100 500 × 500 1000 × 1000
FMM 0.025610 0.001057 0.000266
QFMM (8) 0.025607 0.001057 0.000266
FSM 0.025610 0.001057 0.000266
PDE (non frozen) 0.050516 0.002060 0.000509
PDE (frozen) 0.025610 0.001057 0.000266
PDE (WENO, non frozen) 0.002087 0.000109 0.000023
PDE (WENO, frozen) 0.000004 0.0000001106 0.0000000552
Table 2: Run times in seconds
Method 100 × 100 500 × 500 1000 × 1000 FMM 0.0076630 0.2351160 1.3486670 QFMM (8) 0.0057270 0.1703390 1.0773680 FSM 0.0017270 0.0475360 0.1901890 PDE 0.0895540 11.070000 84.069999 PDE WENO 4.7100000 617.42000 4608.640000
5
Conclusion
As is visible from figures 5 trough 9 the first order algorithms produce the same result both in theory and practice; the difference is negligible. The exception is of course the QFMM with few bins (figure 13) in accordance with the decrease in theoretical accuracy due to heavier quantization during sorting.
2
0 1 2 3 4 5 6 7 8 9 10 x 105 0 0.2 0.4 0.6 0.8 1 1.2 1.4x 10 −6
number of grid cells
time per grid cell (s)
FMM QFMM(8) QFMM(4) FS 0 1 2 3 4 5 6 7 8 9 10 x 105 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x 10 −3
number of grid cells
time per grid cell (s)
First order WENO
Figure 4: Timings of the different distancing methods amortized per pixel as a function of the total number of grid cells.
It is interesting to note the dramatic difference in the PDE based methods when the interface is locked down or not (figures 8, 9, 10, 11). This error source and facile amendment we have not been able to find references to in the literature, but is analogous to the other methods (Q/FMM, FSM) that perform a similar locking procedure. It is clear from the results that the fast sweeping method is the fastest distancing method, in our tests it has the shortest run time, see table 2, and scales optimally (linearly) with increasing grid size, see figure 4. It also has the additional benefit of having the smallest memory footprint as it does not use any auxiliary data-structure or buffers. Second best is the quantized fast marching followed by the normal fast marching. To our surprise the quantized fast marching even though having rid the expensive heap still seems to scale non-linearly (see figure 4). The PDE based versions are notably slower and especially the first order one is not very competitive compared to its fast siblings. When using WENO differentials on the other hand accuracy that is very hard to achieve with the the faster methods can easily be obtained. Figure 12 shows the extreme accuracy, observe that the errors are located at the boundaries and near the singular point in the middle. At the rest of the domain the error is very close to zero.
6
Acknowledgments
This report is the result of an extended project report from a summer course at IDA, Link¨oping University given by professor Ken Museth.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 5: Error of the first order fast marching method. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000.
References
[1] D. Adalsteinsson and J. A. Sethian. A fast level set method for propagating inter-faces. J. Comput. Phys., 118(2):269–277, 1995.
[2] R. Brown. A fast o(1) priority queue implementation for the simulation event set problem. Comm. ACM, 31:1220–1227, 1988.
[3] P. E. Danielsson. Euclidean distance mapping. Comput. Graphics Image Processing, 14:227–248, 1980.
[4] Guillermo Sapiro Liron Yatziv, Alberto Bartesaghi. O(n) implementation of the fast marching algorithm. Journal of computational physics, preprint, 2005.
[5] X.-D. Liu, S. Osher, and T. Chan. Weighted essentially nonoscillatory schemes. J. Comput. Phys., 115:200–212, 1994.
[6] K. Museth, D. Breen, R. Whitaker, S. Mauch, and D. Johnson. Algorithms for interactive editing of level set models. Computer Graphics Forum, accepted for publication, 2005.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 6: Error of the first order quantized fast marching method. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000. 8 discrete bins were used for the quantization.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 7: Error of the first order fast sweeping method. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 8: Error of the first order PDE based distancing method. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000. 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 9: Error of the first order PDE based distancing method. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000. The interface was locked during distancing.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 10: Error of the 3rd order PDE based distancing method. Finite differences are calculated with a HJ WENO scheme, see text. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 11: Error of the 3rd order PDE based distancing method. Finite differences are calculated with a HJ WENO scheme, see text. z shows the error at pixel x,y measured in distance units. Resolution increases from left to right, 100 × 100, 500 × 500, 1000 × 1000. The interface was locked during distancing.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 x 10−3
Figure 12: zoom in on the HJ WENO 500 × 500. Note the magnitude of the error.
0 20 40 60 80 100 0 20 40 60 80 100 0 2 4 6 8 10 x 104 0 20 40 60 80 100 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 4 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 13: Error of the first order quantized fast marching method as the number of bins increase. z shows the error at pixel x,y measured in distance units. The resolution is 500 × 500. The number of bins increase from left to right: 1, 2, 3, and 4 bins respectively were used in the quantization.
0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 x 10−3 0 20 40 60 80 100 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1
Figure 14: The difference between the first order quantized fast marching method and the normal fast marching at 4 (left), 8 (middle), and 1000 bins (right). z shows the difference at pixel x,y measured in distance units. The resolution is 500 × 500.
[7] S. Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer, Berlin, 2002.
[8] S. Osher and J. Sethian. Fronts propagating with curvature-dependent speed: Algo-rithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79:12–49, 1988.
[9] J. A. Sethian. A fast marching level set method for monotonically advancing fronts. Proc. of the National Academy of Sciences of the USA, 93(4):1591–1595, February 1996.
[10] J. A. Sethian. Level Set Methods and Fast Marching Methods. Cambridge University Press, Cambridge, UK, second edition, 1999.
[11] Luminita A. Vese and Tony F. Chan. A multiphase level set framework for im-age segmentation using the mumford and shah model. International Journal of Computer Vision, 50(3):271–293, 2002.
[12] H.K. Zhao. Fast sweeping method for eikonal equations. Mathematics of Compu-tation, 74:603–627, 2004.
