Distance Transforms

<< Click to Display Table of Contents >>

Navigation:  Building Blocks of Spatial Analysis > Distance Operations >

Distance Transforms

In the image processing world a procedure similar to Accumulated Cost Surface (ACS) has been developed in connection with work on pattern recognition. The procedure is known as a Distance Transform (DT), and is typically implemented on binary images in a scanning manner rather than the spreading approach of ACS. Extremely efficient algorithms have been developed to generate exact Euclidean DTs and excellent approximations to exact Euclidean DTs (EDTs) on large images. Some of these procedures have even been implemented within hardware chipsets for real-time applications. GIS vendors do not generally utilize DT methods, but they are included in many image processing toolsets and packages, including the MATLab Image Processing Toolkit: e.g. see the bwdist() function, which provides a variety of metrics for black and white images, including City Block and exact Euclidean. In a series of publications, de Smith (2004, 2006a) has shown how Distance Transforms can be applied to solve complex optimization problems in the geospatial domain. The technique has been applied to a range of practical optimization problems by a number of authors, including selection of optimal routes in road, rail and pipeline engineering projects, some of which are described later in this section.

Using DTs, distances over the image or grid are calculated in an incremental manner based entirely on the distance to nearby (not necessarily adjacent) cells. The standard algorithm involves a two-pass scan of a square or rectangular lattice: a forward scan from top left (NW corner) to bottom right (SE corner), and then a backwards scan from bottom right to top left (i.e. SENW corner). The scan commences by placing the forward scan mask over the first available cell in the grid to be transformed – in a 3x3 scan mask (MxM cells, where M=3) this will be grid cell (row 2, col 2).

A basic DT 3x3 mask is illustrated in Figure 4‑63, where 5 values are used based on the two-value (3,4) integer mask. (3,4) is useful as it is extremely fast to process, and dividing through by 3 at the end yields direct steps of 1 and diagonal steps of 1.33... giving a reasonable approximation to local Euclidean measure (1.414...) and within 6% of exact global Euclidean measure. Typically the underlying grid to be transformed is an array in which unreached cells are initially assigned a large number (e.g. 999) and target or source cells are assigned a value of 0. Cells that act as barriers or restricted areas (e.g. buildings) can be coded with even larger numbers (e.g. 9999), in which case they are never altered. On completion of the two-pass scan each cell in the resulting lattice will contain the distance to the nearest point in the set of source points. If barriers are included they must be in blocks of at least (M‑1)x(M‑1) cells, and accessible areas should consist of at least (M+1)x(M+1) cells. In addition in this case two scan pairs: (NWSE/SENW) and (SWNE/NESW), with one or more iterations, are required to ensure full scanning of the grid. On completion of the scan each cell in the resulting lattice will contain the distance to the nearest point in the set of source points. The pushbroom algorithm described by Eastman (1989) and used in Idrisi for their cost distance modeling is very similar to this approach.

Figure 4‑63 3x3 Distance transformation – scan elements

4

3

4

 

 

 

 

4

3

4

3

0

 

 

0

3

 

3

0

3

 

 

 

4

3

4

 

4

3

4

Forwards

Backwards

 

Combined

The central function in this algorithm is of the form: d0=min{d+D(i),d0}, where d0 is the current value at the central point (0) of the template, D(i) is the local distance to the ith element of the template, and d is the current value at (r,c). The basic algorithm involves O(Mn2) computations where n is the maximum dimension of the grid. Details of DT procedures are provided in the references cited, and code samples in MATLab and Python are available from the author. A selection of 3x3 masks is shown in Table 4‑10.

Table 4‑10 3x3 Chamfer metrics

Case

Description

Adjacency matrix

1

Distances are determined by the L1 or ‘city-block’ metric and paths correspond to the “rook’s move” in chess parlance.

Maximum error: 29.29%

2

1

2

1

0

1

2

1

2

2

Distances are defined as case 1 for any immediate neighbor, i.e. horizontal or diagonal (“rook’s move” and “bishop’s move”).

Maximum error: 41.41%

1

1

1

1

0

1

1

1

1

3

As per 2 but with diagonal distances determined by the Euclidean metric applied locally — sometimes referred to as Octagonal or Local Euclidean metric.

Maximum error: 7.61%

2

1

2

1

0

1

2

1

2

4

Chamfer (3,4)/3 metric. This uses integer values to provide a better estimate of Euclidean distance than Cases 1 or 2.

Maximum error: 6.07%

4

3

4

3

0

3

4

3

4

5

Fractional chamfer (Borgefors, 1986) – non-integer values for both rook’s move and bishop’s move directions (revised after Butt and Maragos, 1998).

Maximum error: 3.96%

1.36039

0.96194

1.36039

0.96194

0

0.96194

1.36039

0.96194

1.36039

These are known as chamfer metrics because the locus of the metric generates a figure similar to a cross-section of a piece of wood with chamfered or beveled edges. The generic form of such masks is defined by the two unique mask values (a,b), e.g. (2,1) or (3,4). The maximum error percentages shown in Table 4‑10 provide the global maximum errors for the metric as compared with global Euclidean measure. The use of 3x3 local metrics is fast and simple, but results in directional bias in the resulting DT. These show the familiar 4- or 8-direction patterns seen in many GIS processes that rely on 3x3 neighborhoods (see, for example, Figure 4‑34).

This problem can be considerably reduced by extending the neighborhood to a 5x5 mask (M=5), in which case 3 values define the mask: (a,b,c) as shown in Figure 4‑64. Amongst the best integer values for this mask size are the triple (5,7,11), dividing all cell values by 5 at the end to give (1,1.4,2.2). The best fractional values are (0.9866, 1.4141, 2.2062). Use of the integer values above ensures that every grid cell distance computed (i.e. global values) is within 2% of the exact Euclidean distance and within 1.36% for the fractional masks. This compares with global errors of up to 7.6% with a pure local Euclidean metric.

Figure 4‑64 5x5 DT mask

2b

c

2a

c

2b

c

b

a

b

c

2a

a

0

a

2a

c

b

a

b

c

2b

c

2a

c

2b

Figure 4‑65A provides the simplest test case for the DT procedure. This shows a uniform region with a single target point in the center (originally coded as 0, with all other grid cells coded 999). This figure shows the best non-exact DT possible using a fractional-valued 5x5 scanning template. The exact Euclidean distance transform (EDT) of this uniform space yields a perfectly circular set of filled contours around the central point. The vector equivalent would be a multi-ring buffer operation. When applied to multiple points (e.g. as part of a point-pattern analysis exercise) the set of computed distances, F(), is known as the empty space function.

If a simple barrier is introduced, e.g. representing a river with a bridge point, shortest paths appear diffracted through the gap (Figure 4‑65B). In this example the target point has been set as (2,50). With a complex barrier, as in Figure 4‑65C, a single extra iteration has been used to ensure full convergence. In all cases shortest paths to target points can be drawn as orthogonal to the distance bands, but are best constructed using tracking of the selected steps, as described further in the examples below. Although the DT approach may seem somewhat divorced from practical GIS analysis, DT procedures have been extended and developed by: Ratti et al. (architecture and the built environment; 2001, 2004); de Smith (transport routes; 2004, 2006), Ikonen and Toivanen (image processing; 2005); Kristinsson et al. (geothermal pipelines, 2005); and Wei Li et al. (railways in mountainous regions, 2017) to provide a new and fast family of ACS-like procedures for solving many complex grid-related distance problems.

Figure 4‑65A Distance transform, single point. 5x5 DT, no constraints, target (50,50)

clip0116.zoom81

Figure 4‑65B Distance transform, single point 5x5 DT, barrier – target (2,50)

clip0117.zoom78

Figure 4‑65C Distance transform, single point. 5x5 DT, complex obstacle – target (50,90)

clip0118.zoom67

We limit discussion here to a small sample of applications. The first deals with radially symmetric models of traffic in cities ― see Angel and Hyman, Urban Fields (1976, pp 15-16 and 159) for a full discussion of the models involved. The second example examines pedestrian movement within urban areas, whilst the third set of examples address questions of optimal road, rail and pipeline location under a range of constraints.

Figure 4‑66 shows the result of applying a least ‘cost’ distance transform (LCDT) to a city with traffic speeds symmetrically varying from almost 0kph in the center (where the crossed lines meet) to much higher levels (e.g. 50kph+) at the city edge. The cost field in this case is traffic speed (i.e. a velocity field) and the curved lines show contours of equal travel time (isochrones) to or from a point due South of the city center.

Figure 4‑66 Urban traffic modeling

clip0119.zoom70

The central function in the DT procedure in this instance is of the form:

d0 = min{d+D(i)*COST(r,c), d0}

where d0 is the current value at the central point (0) of the mask, D(i) is the local distance to the ith element of the mask, d is the current value at (r,c) and COST(r,c) is the (averaged) cost function at row r, column c. As can be seen, the velocity field distorts the distance bands. Shortest (quickest) paths across this transformed surface may be plotted by following routes to or from the source point in any given direction, remaining at right angles (orthogonal) to the equal time contours, or by tracking of the computed quickest paths. Hence the quickest route from the start point due South of the center to locations North of the city center would curve around (East or West) the center to avoid congestion (low traffic speeds. In this particular example these curve form log spirals (Wardrop, 1969). Angel and Hyman applied models of this type to the analysis of traffic flows in Manchester, UK. Distance transform methods enable such models to be calibrated against known analytical results and then applied to more realistic complex velocity and cost fields, for example to examine the possible effects of introducing central zone congestion charging.

A second example of this procedure is shown in Figure 4‑67. This shows a section of the Notting Hill area of central London, with the annual Carnival route shown in white and areas in dark gray being masked out (buildings etc.). Colors indicate increasing distances (in meters) from the Carnival route along nearby streets.

Figure 4‑67 Notting Hill Carnival routes

clip0120

Although this example has been created by a DT scanning procedure the result is sometimes described as a “burning fire” analogy, with the fire spreading from the source locations (in this case a route) outwards along adjacent streets and squares (for true fire spread models see www.fire.org and Stratton, 2004). The same techniques can be applied in public open spaces and even within building complexes such as shopping malls, galleries and airports. This illustrates the ability of DT methods to operate in what might be described as complex friction environments.

Distance transforms are very flexible in their ability to incorporate problem-specific constraints. For example, in a typical ACS procedure variable slopes are modeled by regarding steeper slopes as higher cost, rather than relating the spreading process to slope directly (other than some uphill/downhill constrained implementations). The selection of a route for a road, railway or pipeline is affected by slope, but primarily in the direction of the path rather than the maximum slope of the physical surface (although this does have a bearing on construction cost).

Figure 4‑68 illustrates the application of a gradient-constrained DT (GCDT) to the selection of road routes through a hilly landscape. The most northerly two routes (upper right in the diagram) were constrained to a 1:12 and 1:11 gradient (before engineering grading) whilst the most southerly two represent routes found using less restrictive 1:9 and 1:8 gradient constraints.

Figure 4‑68 Alternative routes selected by gradient constrained DT

clip0121.zoom69

This technique has been applied to the problem of mountain railway alignment optimization (Wei Li, Hao Pu, Paul Schonfeld et. al, 2017). As the authors explain

"Various complex constraints must be satisfied in railway alignment design for topographically complex mountainous regions, where the natural terrain gradient between the start and end points greatly exceeds the maximum allowed design gradient. The alignment design for such environments is so challenging that existing methodologies have great difficulties in automatically generating viable railway alignment alternatives. We solve this problem with an alignment optimization methodology combining the bidirectional distance transform and a genetic algorithm. This methodology consists of a two-stage search that first identifies promising corridors within a broad area and then optimizes alignment alternatives within those corridors.... This approach has been employed by the China Railway Eryuan Engineering Group and has been applied in designing the Sichuan-Tibet, Tongdao-Quanzhou, Zhengzhou-Guiyang, and Ganluo-Yaan railways as well as a railway across South America in Peru and Brazil. "

In the case study area illustrated below (Figure 4-69A, B), an area of extremely mountainous terrain covering over 900 sq km is shown. This forms part of the Sichuan-Tibet railway design project. The blue line shows the Distance Transform optimized path, which differs markedly from the manually aligned path chosen by specialists in the field, and which surprisingly crosses itself. The reason for this unexpected behavior is that the constrained path does not lie on the terrain surface, but almost entirely under the surface (in tunnels) or above the terrain surface via bridges. Where the path crosses itself there is a large difference in the vertical location of the relevant sections of the path, which occur because of the gradient constraints applied. The final optimized path, generated using genetic algorithm (GA) techniques with the DT path as a starting point, is shown in Figure 4-69b. It represents an estimated saving of approximately $3million per annum over the best manually aligned path.

Figure 4-69A Sichuan-Tibet railway design project — initial route selection

Fig.11 Optimized and close-to-manual alignment paths generated by DT

Figure 4-69B Sichuan-Tibet railway design project — optimized route

Fig.15-a Horizontal alignment

With minor modifications the GCDT procedure described above has been utilized with a high resolution (10meter) DEM to locate alternative pipeline routes for the Hellisheiði 80MW geothermal power station in Iceland (Kristinsson, Jonsson and Jonsdottir, 2005). In this case additional constraints were applied giving the maximum and minimum heights permitted for the route and differential directional constraints that permitted limited flow uphill, which is possible for these kinds of two-phase pipelines (Figure 4‑70A,B).

Figure 4‑70 Hellisheiði power plant pipeline route selection

A. Overview of power plant area

clip0122

B. Alternative solution paths A, B1, B2 and C for Well 3

clip0124