Geodesic distances

The geodesic distance between two pixels p and q is defined as the length of the shortest path from p to q. Suppose $P=\{p_1,p_2,...,p_n\}$ is a path in a connected domain between pixels p₁ and p_n, i.e. p_i and p_i+1 are connected neighbors for $i \in \{1,2,...,n-1\}$ and p_i belong to the domain for all i. The path length l(P) is defined as

$$l(P) = \sum_{i=1}^{n-1}d_N(p_i,p_{i+1})$$ (2.30)

the sum of the neighbor distances d_N between adjacent points in the path. For instance, the city-block metric dist_CB defined on the rectangular image domain is a geodesic distance where only 4-direct neighbors are connected and d_N=1. Similarly, the chess-board metric dist_chess considers that 8-direct neighbors are connected and d_N=1. The chamfer metric dist_cha(3:4) considers 8-direct neighbors to be connected, d_N=3 between 4-direct neighbors, and d_N=4 between diagonal neighbors. Finally, the chamfer 5-7-11 considers that all neighbors within a $5 \times 5$ neighborhood are connected, with $d_N=5,7 \; \mathrm{or} \; 11$. On the other hand, with the Euclidean metric dist_e, there is a direct path between any two pixels in the image, which makes the problem non-local and renders it so much more difficult to compute.

In the previous sections, distances are always defined for every pixel of the image, i.e. on a rectangular domain. In contrast, Piper and Granum [122] study geodesic distances defined on any connected domain, both convex and non-convex ones. On convex domains, they show that the usual two raster scans algorithm of section does not function anymore. As illustrated at Figure , the chamfer masks do not support propagation in every direction, but only in two 135^o angles.

  

Figure 2.13: Left: chamfer masks and the supported propagation directions. Right: A convex domain on which the two raster-scan chamfer DT algorithm does not compute the distance transformation

The surprising part should actually be that this algorithm does work on rectangular image domains. The reason for this is that, in a geodesic DT, one only needs to support the propagation along one of the paths of minimal length to get the correct distance value. For any pixel located between the up-right and right directions from their nearest object pixel, the two raster scan chamfer algorithm does provide one (and only one) such path, made of right steps only, followed by up-right step only. For arbitrary convex domains such as that of Figure , this particular path does not always belong to the domain, and another minimal-length path should be used instead, which is impossible given the available supported propagation directions.

There are several solutions to this problem. First of all, one could apply the two-scan algorithm iteratively until no pixel changes value anymore. Unfortunately, there is no way to predict how many iterations should be used. Instead, they propose to use 4 different raster-scanning orders with either 3 or 4 neighbors in each mask, or only 3 different raster-scans and masks similar to those used by Ragnelmam [128] at Figure .

For non-convex domains (Figure ), even 4 raster scans algorithms do not always reach the whole domain, and would need to be applied iteratively. Instead, Piper and Granum suggest to use propagation algorithms, starting from the object pixels, then considering their neighbors, their neighbors' neighbors, ... They propose two simple implementations of such algorithms, recursive propagation which is depth-first, and ordered propagation which is breadth-first. The ordered propagation algorithm they propose uses the neighboring order, i.e. the geodesic chess-board distance order. This is suboptimal when another metric is used.

  

Figure 2.14: Left: a non-convex domain. Right: distance transformation using the geodesic equivalent of the dist_chess metric.

Verwer, Verbeek and and Dekker [167] propose an efficient algorithm for ordered propagation. The main idea is to scan the pixels in the domain in the order defined by the metric used for the DT. With a geodesic metric, where the minimal path length is defined as equation () with d_N > 0, this guarantees that every pixel will only be propagated once. Therefore the computational complexity of the DT is optimal, i.e. strictly proportional to the number of pixels.

Scanning the pixels in the metric order is made possible by sorting the pixels in the propagation front before propagating them. With the metrics we consider, there is only a finite number of possible distance values. Then, bucket sorting has an optimal O(n) computational complexity where n is the number of elements to sort.

Practically, there is one bucket bucket(d) per possible distance value d. Initially, every object pixel is put into bucket(0), all other buckets are empty. D(p) is set to "0" for every object pixel and to "$\infty$" for every non-object pixel in the domain. In the propagation phase, illustrated at Figure , pixels are taken from the bucket with the smallest value d for which bucket(d) is not empty. For every neighbor q of the current pixel p, we suppose that D(q)=D(p)+d_N(p,q), i.e. that the minimal-length path (eq ) from the object to q goes through p. If this leads to a smaller D(q) than currently stored, its value is updated and q is inserted in bucket(D(q)). The algorithm stops when all buckets are empty.

  

Figure 2.15: Bucket sorting propagation.

The implementation of the above algorithm requires special care. Because one cannot know a priori the size of the buckets, the memory for those has to be dynamically allocated. If this is done with classical dynamic lists, with one pixel per element, there is at least one memory allocation operation per pixel in the image. Dynamic memory allocation then becomes the major factor in the computational cost. Instead, memory should be allocated by chunks, with a number c of pixels stored in each chunk. Verwer analyzes the influence of the chunk size c on the memory requirements and processing time. He concludes that any chunk size between 15 and 100 is satisfactory. With c<15, too much processing time is wasted in memory allocation. With c>100, too much memory is wasted in partially filled chunks. This upper limit is of course dependent on the image size.

Also, during the propagation, all buckets are never used simultaneously. Actually, when processing bucket(d), only buckets between d and d+max(d_N) can contain pixels. Therefore, Verwer proposes to reuse the buckets circularly.

Provided it is carefully implemented, Verwer's algorithm appears to be optimal for geodesic distances defined with path lengths as in ().