{'fraction_non_alphanumeric': 0.036761898149755265, 'fraction_numerical': 0.01829416461276773, 'mean_word_length': 4.168281306063868, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Distance Transform and its Computation

Distance transformation is an image processing technique used for many different applications. Related to a binary image, the general idea is to determine the distance of all background points to the nearest object point (or vice versa). In this tutorial, different approaches are explained in detail and compared using examples. Corresponding source code is provided to facilitate own investigations. A particular objective of this tutorial is to clarify the difference between arbitrary distance transforms and exact Euclidean distance transformations.

1 General description 1.1 Examples and possible applications Figure 1 shows what can be achieved in general by a distance transform. The input is typically a binary image (Figure 1a), i. e. the pixels can have only one out of two different values. One value is associated to the background and the other value defines the object pixels. The object pixels can belong to one or more (disconnected) objects. Depending on the application, either the distances inside the objects are of interest (Figure 1b) or the distances outside the objects (Figure 1c).
The segmentation of images in regions is one important application of the distance transform. Afterwards, each region represents one object point or a cluster of connected points. This is highly related to so-called Voronoi diagram. Figure 2 and Figure 3 shows two examples. From these Voronoi diagrams it can be deduced which points are adjacent. This information could be helpful in solving graph-based problems.
The distance within objects also could be of interest for finding the centre of objects like hands or for the skeletonization of objects, see Figure 4. Intermediate processing results may vary when different distance metrics (which are discussed in Subsection 1.3) are used, see Figure 5. And last but not least, distance transforms can support the segmentation of overlapping objects like cells, Figure 6. This paper wont discuss, how the distance transform has to be used in achieving the goals of these applications; instead, it concentrates on the basics of distances and on algorithms for computing a distance transformation in an appropriate manner. First ideas about distances in digital images have been discussed in a pioneering work by Rosenfeld and Pfaltz in 1968 [Ros68].
The following two Subsections discuss neighbourhood relations and different metrics for distances. If you are already familiar with it, just skip to Section 2.

Neighbourhood relations
In an orthogonal grid of pixels, we can differentiate between two neighbourhood relations. The first is called 4-neighbourhood since only the horizontal and vertical neighbouring pixels are considered. When also the diagonal pixels are taken into account it is called 8-neighbourhood, This differentiation has implications in several directions. If, for example, binary image objects have to be counted, the result can be different. Figure 7 b) shows dark object pixels, which belong to a single object. In 8-neighbourhood, all pixels belong to a single object, because the marked pixels are considered as neighbours. However, in 4-neighbourhood the marked pixels are not connected and there are two binary objects.

Minkowski metric
The distances shown in Figure 1 are based on the Euclidean metric. This is what we also know as beeline distance. If we hammer two nails into a board and connect them with a rubber band, then the Euclidean distance is equal to the one-way length of the rubber band. In two dimensions this distance is defined by the Pythagorean theorem: c 2 = a 2 + b 2 . If there are two points p and q with two-dimensional coordinates (p x , p y ) and (q x , q y ), then the Euclidean distance d is computed with d = (p x − q x ) 2 + (p y − q y ) 2 , see Figure 8. Depending on the application, different metrics could be useful. A more general definition is given by the Minkowski distance: with D being the number of dimensions. The Euclidean metric is a special case with e = 2. Alternatively, e = 1 corresponds to the city-block distance, see Figure 9. This distance is sometimes called the Manhattan distance because most streets in Manhattan are either parallel or orthogonal to each other. If somebody wants to walk from point p to point q, then she has to go along the streets and there is no possibility to cross the buildings.
Setting variable e to infinity defines the chessboard distance. This is the number of moves that the king has to take on the chessboard from one square to another, Figure 10.

Chamfer metric
Apart form metrics that can be derived from the Minkowski distance in (1), there are other possibilities to define distances. One is called 'chamfer' distance and it is constructed by using a different step size in diagonal directions compared to horizontal or vertical direction. Both, the city-block and the chessboard distances apply steps of '1' in horizontal and vertical direction. However, the diagonal step is equal to '2' in the former and equal to '1' in the latter case as can be seen in Figures 9 b) and 10 b). As the Euclidean distance would fit many applications best, there have been made lots of attempts to approximate this metric by simpler constructions without the need of multiplications and drawing the square root. So, one of the first ideas was to use a diagonal step of √ 2. The resulting effect can be seen in Figure 11 a). The distance values directly close to the centre pixel, in horizontal, vertical, and diagonal direction are now exactly the same as for the Euclidean distance (Figure 8b). At other positions, the distances differ. For example, the distance measure at the pixel in the bottom right corner deviates from the Euclidean distance by ∆d = √ 44.3 − √ 41. Borgefors has derived in [Bor84] that the average and the maximum approximation error can be minimized by using a diagonal step of 1.351. That means, the optimal relation between diagonal and other steps is Using the approximation 4/3, the distance computations can now be performed in integer arithmetic if all values are scaled by factor 3. Figure 11 b)

Non-Euclidean Distance Transforms
This section explains several versions of distance transformations that can be implemented by fast algorithms but are not able to produce exact Euclidean distances. If you are only interested in exact Euclidean distance transformations, then you can simply skip to the section 3. In the following, the case shown in Figure 1 c) is considered. That is, the shortest distances between a background pixel and any object pixel are searched for.

City-block distance transform
One of the earliest ideas to perform a distance transformation was based on distance propagation. Starting from the positions of the object pixels, whose distances are zero by definition, the distances of the background pixels are determined step by step while incrementing the distance value. This can be achieved by sequentially processing all image rows first from top-left to bottom-right and then in backward direction.
At first, the distance image (or matrix) has to be initialized with values of zero at the positions of object pixels and with sufficient large value (larger than the maximum possible distance) at all other positions. An example is given in Figure 12 a). Then all rows are scanned from left to right. If the distance d(i, j) value at the current position is larger than the predecessor plus one d(i, j − 1) + 1, it is set to a new value d(i, j) := d(i, j − 1) + 1. Afterwards, a) This procedure propagates the distances to the right and to the bottom as shown in Figure 12 b). This direction of propagation has to be reversed in a second scan from the bottomright up to the top-left position: and the final result in Figure 12 c) is obtained. The corresponding source code is shown in Listing 2 (see appendix). In total, the image is scanned twice and two comparisons per scan have to be made at each pixel position. An alternative approach is based on four scans, but with only one comparison per pixel access. It first processes all columns downward and upward: Since the comparisons are independent of each other, all columns could be processed in parallel. Afterwards, all rows are processed in the same manner:  Figure 14 a) depicts the distance transform result for a larger example with ten isolated object pixels. The brighter a pixel is in this picture, the higher is its distance to the closest object pixel. The typical diamond shape of the city-block distance propagation can be seen clearly (compare also Figure 9 5 4 3 2 1 2 3 2 3 4 4 3 2 1 0 1 2 1 2 3 3 2 3 2 1 2 1 0 1 2 2 1 2 3 2 2 1 0 1 2 1 0 1 2 2 3 2 1 2 3 2 1 2 2 1 0 1 1 3 3 3 2 3 3 2 1 1 0 1 2 4 3 4 4 3 2 2 1 2 3 5 4 5 5 4 3 3 2 3 4

Approximate Euclidean Distance Transform
Danielsson had the idea to not propagate the distances but the relative vertical and horizontal positions [Dan80]. Object pixels always have a relative position of [∆ y , ∆ x ] = [0, 0], see Figure 16 a). Each horizontal step increments the horizontal relative position ∆ x and each vertical step increments the vertical relative position ∆ y . A diagonal step increments both. The Euclidean distance can be computed with d = (∆ 2 y + ∆ 2 x ) 0.5 . The source code in Listing 4, Listing 5, and Listing 6 differs somewhat from Danielsson's proposal and follows the modifications of [Gre07].
In addition to the distance matrix distMat that has to be generated, two other matrices distMatY and distMatX have to be maintained. They store the propagated relative positions. Based on these relative positions, the squared Euclidean distances have to be determined. Instead of computing these distances again and again, a function getDist() is used. It calculates the different distances only once, stores these values in a matrix preCalc, and also utilizes the symmetric property of the metric.
The source code in Listings 5 and 6 is somewhat lengthy, since several decisions have to be made. The processing comprises two scans: the first is starting in the second row going down to the bottom row and the second is operating in opposite direction. In each scan, the rows are processed forward and backward. The lines 3-14 in Listing 5, for example, perform a vertical propagation of relative positions. The values from the top neighbour position are taken (yr = distMatY(i-1,j); xr = distMatX(i-1,j);). The vertical component is incremented (yr+1) and the corresponding squared distance is determined using the function getDist(). If the current distance is longer than this new one (if distMat(i,j) > dist), it is replaced and also the relative positions are updated (distMat(i,j) = dist; distMatY(i,j) = yr+1; distMatY(i,j) = xr). However, function getDist() may not be accessed when the relative positions are larger than the dimensions of matrix preCalc (if yr < maxDist).
This procedure is repeated for the two other positions in the causal neighbourhood (i,j-1) and (i-1,j-1). The subsequent backward scan for this row checks the positions (i,j+1) and (i-1,j+1). Figure 16 b) shows the result of this first scan. Note that the sorted pairs a, b and b, a correspond to the same distance d 2 = a 2 + b 2 . However, the derived order of a and b indicates to which object pixel each background pixel has been assigned.
The second scan propagates the minimum distances from the bottom row upwards (Listing 6) in a very similar manner and produces the final result (Figure 16 c). For the given example, this corresponds to the exact Euclidean distances as shown in Figure 18.
Under special circumstances, it is impossible to derive the correct distance. This is because the algorithm depends on the propagation of distances, while the assignment of the closest object pixel can lead to disconnected background pixels. This has already been mentioned and explained by Danielsson himself [Dan80] and different numeric examples have been given in [Bai05] and [Fab08]. One example can be seen in Figure 17. The continuous space is separated in regions by drawing the perpendicular bisector between two points. There is a clear assignment of each position to one of the three point, Figure 17 a). In discrete space, the three object pixels (marked with a cross) also split the set of background pixels in three region, Figure 17 b). The white lines show the separation in continuous space. The discretization causes a disruption of the middle region. The approach of Danielsson is not able to take account of this disruption and assigns the isolated pixel to the bottom region with a relative position of 1, 13. There is no chance to propagate the positions from '4,10' to the correct relative distance of '5,12' via '5,11' or via '4,11' because the corresponding positions are closer to the other object pixels.
A typical visual appearance of this distance transform result is shown in Figure 14  The bright isolated pixel has been assigned to the wrong object pixel with a distance of d 2 = 1 2 + 13 2 = 170. However, its correct distance is d 2 = 5 2 + 12 2 = 169 with relation to the object pixel in the middle.

Exact Euclidean Distance Transform
Some applications depend on the exact calculation of the distances between each background pixel and the nearest object pixel. Thus, the following explanations refer to the case shown in Figure 1 c). In principle, one could first identify the inner contours of all objects. The inner contour includes all object pixels that have at least one background pixel in their neighbourhood. In a second step, one could compute the Euclidean distances from each background pixel to all contour pixels and find the minimum. This minimum distance is then stored at the position of the examined background pixel. In worst case, each object pixel could be a contour pixel 1 . Let n O denote the number of contour pixels. The total number of pixels is given by the image size L × M. The number of background pixels is then n B = L · M − n O . The number n c of distance calculations would be The maximum of required calculations max(n c ) = (L·M) 2 /4 is reached for n B = n O = (L·M)/2. Consequently, this naive approach has a quadratic complexity with respect to the image size: For a better understanding of the following text, the most important variables are summarized in Table 1.

Basic approach
Many researchers have thought of algorithms that can compute the exact Euclidean distance transform (EEDT) in a faster manner. Some of these successful methods are based on the separate processing of all dimensions [Mei00, Mau03,Fel12]. Let us take an example with six object pixels. The distance matrix in Figure 18 a) shows the positions of object pixels with a value of zero. The following derivations use squared Euclidean distances . . index of current or particular row i n . . . vertical coordinate of an object pixels p n j . . . column index j 0 . . . index of current or particular column j n . . . horizontal coordinate of an object pixels p n j s . . . horizontal location of an intersection, j s ∈ IR k . . . column index related to object pixels L . . . number of columns M . . . number of rows n . . . index for object pixels p n . . . object pixel a) with p = (p x , p y ) being the position of a background pixel and q = (q x , q y ) being the position of an object pixel. The square root can be drawn (if necessary) after all squared distances have been determined. The distances for the background pixels are not known yet and they are initialised with a sufficiently large value higher than (M − 1) 2 + (L − 1) 2 , which represents the maximum possible distance for a M × L matrix. In Figure 18, this value is indicated by '∞'.
From (4) it can be seen that the distance consists of two additive terms, which can be computed independently. For example, one could first calculate all vertical components D v = (p y − q y ) 2 . The result is depicted in Figure 18 b).
In a second step, the horizontal components D h = (p x − q x ) 2 have to be determined row by row. At all positions in each row, it has to be checked which combination D = D v + D h of vertical and horizontal components results in the shortest distance. Starting in the topleft corner, there is no distance assigned yet. The row is scanned to the right for the next available distance component. A value of D v = 16 can be found at a horizontal distance of d h = 1. The combination of both components is D 2 = 1 2 + 16 = 17. The index k in D k indicates the column of the used vertical distance. This is the first candidate. The second can be computed at a horizontal distance of 4: D 5 = 4 2 + 1 = 17. Two more candidates are available: D 6 = 5 2 + 25 = 50 and D 8 = 7 2 + 4 = 53. The best candidate is obviously This value is written to the current pixel position. The candidates for the second pixel in the first row are: D 2 = 0 2 + 16 = 16 (pure vertical distance) D 5 = 3 2 + 1 = 10 D 6 = 4 2 + 25 = 41 D 8 = 6 2 + 4 = 40 .
The already existing value of 16 has to be replaced by D * = 10, because this is the shortest distance to any object pixel. This procedure continues until all distances are fixed for this row. Then all other rows are processed in the same manner.
Listing 7 contains a possible algorithm for this procedure. This source code allows an estimation of the corresponding complexity. Lines 15 to 34 compute the vertical components. Interestingly, there is no multiplication involved despite the fact that (p y − q y ) · (p y − q y ) has to be computed for each background pixel. The outer loop for j=1:L steps through all columns. There are two inner loops: for i=2:M and for i=M-1:-1:1. The first one scans the current column downwards starting with the second row and the second loop scans the same column upwards. Both try to propagate distances in such a manner that larger distance entries are overwritten by possibly smaller ones.
With respect to the example in Figure 18 a), the downward scan firstly enfolds its impact in the second column (counted from the left) at the position right under the object pixel with distance equal to zero. At this position it finds a distance value that is larger (= ∞) than the distance above plus the actual distance step, which had been initialized with distStep=1. The current value is set to 0 + 1 = 1 and the distance step is increased by two. What's the deal with that increment?
Given two consecutive numbers x and x + 1, the difference between their squared values is (x + 1) 2 − x 2 = 2x + 1. Starting with x = 0, the first difference is equal to one (distStep=1). Then the difference between the squared values grows by 2 (distStep = distStep + 2;) with each increment of x. Using this step as an additive term in distMat(i,j) = distMat(i-1,j) + distStep, the algorithm computes the required squared distance values without any multiplication.
As soon as the position of an object pixel is reached, the current distance value is not larger than the propagated one; the distance value remains unchanged and distStep is reset to one. The result of the downwards propagation is shown in Figure 19. The upward scan within a column does not change anything before the first object pixel is passed. After this position the distances are update if they are larger than the propagated values, see 18 b). These scans are rather fast and there is only opportunity for improvements in the upward scan, if there are many columns without any object pixel (entire column could be skipped in the upward scan) or The current minimum distance is D * related to point p 2 from column k * . The distance measure d 3 cannot be shorter because its horizontal component ∆ x is already longer than the current minimum: ∆ x > d * .
the object pixels concentrate in the upper rows of the image (lower part of the column could be skipped). In maximum there are L · M distances to be computed. So, the vertical propagation has linear complexity with respect to the image size. The processing of rows needs more efforts and defines the complexity of the entire approach, Listing 7 lines 37-52. The outer loop for i=1:M takes care of all matrix rows. The corresponding row of vertical distances is copied to a new vector distMatV because the values in matrix distMat can be overwritten. In lines 40 and 43, there are two nested inner loops, both scanning all positions within the current row. Index j of the first inner loop defines the position for which the minimum distance has to be found. The second inner loop scans for all possible candidates of vertical distances distMatV(k). Figure 20 shows a situation where a background pixel at position (i 0 , j 0 ) has to be assigned the distance to the closest object pixel p k . Scanning the line i 0 from k = 1 towards k = L, the current point has already been compared to the points p 1 and p 2 resulting to a minimum distance of d * = √ 4 2 + 3 2 = 5 in this example. The next vertical distance D v < ∞ can be found in the column of point p 3 .
In the worst case, all positions j 0 have to be compared with all positions k with 1 ≤ j 0 , k ≤ L, i.e. L·(L−1) distances have to be computed. Taking all rows into account, we have M ·L·(L−1) distances to calculate leading to a O(ML 2 ) complexity, which is at least somewhat lower than for the naive approach mentioned in the beginning of Section 3.

Faster approaches
The most important question is, whether it is possible to skip some of the distance calculations. While the determination of vertical distances is already quite fast and has only linear O(M) complexity, the calculation of distances along the rows has some room for optimizations. With respect to the basic approach discussed above, there are some tricks to speed up the processing.

Change order of column and row processing
Depending on the aspect ratio of the given image, the complexity O(LM 2 ) is lower than O(ML 2 ) if M < L. If the complexity of the row-wise processing will not reduced by another mechanism, the order of processing columns and rows should be exchanged in this case.

Avoid unnecessary computations
The computations in lines 43-49 of Listing 7 should only be executed if there is a chance to find a new minimum distance. Columns without any object point contain infinite (=maxDist) vertical distances. So, checking 'if distMatV(k) < maxDist', it can be ensured that variable k is not pointing at such column and unnecessary computations can be avoided. However, if most of the columns are occupied by at least one object pixel, this has little effect.
A second simple trick that can be combined with the first one concerns the horizontal component of the distances to be computed.
As soon as a position k = j 0 + ∆ x with ∆ x > d * is reached, no shorter path than d * can be found and thus the scan can be terminated. This situation is shown in Figure 20. The current minimum distance d * refers to point p 2 and the next available point p 3 has a horizontal distance that already is larger than d * and the processing of loop for k=1:L can be stopped. This leads to large saving effects when j 0 is small. However, the effect decreases the more the current position j 0 moves towards L. If the scanning process is split into a forward scan j 0 < k ≤ L, which can be processed as described and a backward scan from k = j 0 − 1 to k = 1, then it is possible to stop both scanning processes as soon as |k −j 0 | ≥ d * is reached. These modifications are implemented in Listing 8. Remember that all distances are squared ones.

Fast exact Euclidean distance transforms
The previous Subsections explained how the exact Euclidean distance transform can be realised and they proposed to skip unnecessary computations. However, the discussed methods for doing so are quite simple and have impact only under special conditions. As mentioned in The explanations in Subsection 3.1 regarding the combination of the vertical and the horizontal component to find the final distance in (4) can be viewed from another perspective. Each object point p n located at (i n , j n ) contributes to row i 0 with a vertical component   Figure 21.
Since we are looking for the minimum distance, the minimum value is determined for each column j. This results to a row of distances which is equal to what we already got in Figure  18 c), in first row.
Additionally, we can keep track of which column k contributes the minimum distance for a certain column j 0 . As can be seen in Figure 21, the background pixel at position j = 1 has a minimum distance to the object point lying in column k = 2 and the same distance to the point lying in column k = 5. The background points in columns 2 ≤ j ≤ 6 are closest to the object point from column k = 5. For j = 7, again two object points share by chance the same distance to the background point a.s.o. The object pixel from column k = 6 is not taken into account for this row because already its vertical component D 6,v = 25 is to high.

Method by Meijster, Roerdink, and Hesselink
Looking at Figure 21, it becomes clear that the values in the rows for D k follow a quadratic equation with j being the independent variable. Hence, the distance D k = f (j) can be visualized as parabolas, see idea has been proposed by Meijster et al [Mei00]. Working source code has been provided by Felsenszwalb and Hutterlocher [Fel12] while it is unclear whether they had been aware of Meijsters work. So, the key point is to determine the position of the intersection of two parabolas. Let us assume that there are two parabolas with vertices at the columns k and l. Since they have a single intersection, the position j s of this intersection can be determined using (5) with following equation: Using the example in Figure 22, an algorithm is explained in the following, which can significantly reduce the complexity. Two vectors are needed for storing the required information of the lower envelop. The first vector ks will contain the positions k ∈ {j n } of all parabolas which finally contribute to the lower envelop and the second vector js will contain the positions j from which a parabola starts to contribute to the lower envelop. Theses positions can be derived from the intersections j s rounded to the next larger integer value.
The vectors have to be initialized with the parameters of the first, possibly virtual, parabola: ks(1) = 1 and js(1) = −∞. The huge negative value serves as stopping criterion for the algorithm as will be explained below. The corresponding source code is shown in Listing 9 in lines 6 -12.
The position of the first available (blue) parabola is k = 2 and the intersection j s with the virtual parabola is per definition larger than −∞. The information of the new parabola is stored: ks(2) = 2 and js(2) = max(1, j s ). It has to be ensured that the lower border starts not before j = 1.
The position of the next available (red) parabola is k = 5 and the intersection with the previous is exactly at position j s = 1. Since j s <= js(2) holds, the blue parabola is hidden by the new one and does not contribute anymore. Its parameters can be overwritten: ks(2) = 5, js(2) = 1. This check has now to be done with all previous parabolas (lines 21-26). As soon as the intersection is j s > js(idx), the previous parabola contributes to the envelop and the information of the current parabola is appended to the vectors ks and js. The third (magenta) parabola of the example in Figure 22 intersect with the previously stored (red) parabola at a position beyond L and thus can be ignored. The fourth and last (green) parabola has an intersection with the red one at j s = 7. This is larger than the previously stored value of js(2) = 1, so it does not hide the previous parabola and the vectors must be complemented with the parameters of the green parabola: ks(3) = 8, js(3) = 7 (lines 27 -31). There is no further parabola available and the intersections vector can be closed with js(1) = 1 and js(idx + 1) = L + 1 defining the maximum range. The vectors now read as js = (1, 1, 7, 11) and ks = (1, 5, 8) giving the information that parabola at k = 1 is either not existing or not contributing (1 ≤ j < 1 !), parabola k = 5 defines the distances for the positions 1 ≤ j < 7 and parabola k = 8 is responsible for the distances at 7 ≤ j < 11. Using equation (5), the squared Euclidean distance can be calculated for all positions j (lines 37-42). This procedure is performed for each row of the distance matrix.
In total, each row is scanned twice. The first scan identifies the contributing object points (parabolas) and the second scan applies the collected information and calculates the final distances D(i : 0, j 0 ). Even though, the first scan uses a double while-loop, its complexity can be considered as being linear (O(L)). That means, the entire method has O(M · L) complexity.