IN
DEGREE PROJECT MATHEMATICS, SECOND CYCLE, 30 CREDITS
,
STOCKHOLM SWEDEN 2018
Fast marching and fast
sweeping in optimal path
planning
HÅKON JARVIS WESTERGÅRD
Fast marching and fast
sweeping in optimal path
planning
HÅKON JARVIS WESTERGÅRD
Degree Projects in Scientific Computing (30 ECTS credits)
Degree Programme in Applied and Computational Mathematics (120 credits) KTH Royal Institute of Technology year 2018
Supervisors at NTNU Trondheim: Stein Krogstad, Knut-Andreas Lie Supervisor at KTH: Sara Zahedi
TRITA-SCI-GRU 2018:298 MAT-E 2018:68
Royal Institute of Technology
School of Engineering Sciences
KTH SCI
i
Preface
This Master’s thesis has been conducted at the Royal Institute of Technology (KTH) in collabora-tion with the Norwegian University of Science and Technology (NTNU) in Trondheim in spring 2018. The thesis is part of the program "Master’s Programme in Applied Mathematics". The topic "Level Seth Methods" applied to orienteering map was suggested by the supervisors at SINTEF as the author and one of the supervisors have background as a orienteering runner.
The readers of this project are assumed to be familiar to some basic graph theory, partial differential equations and numerical treatment of partial differential equations.
ii
Acknowledgment
I would like to thank Stein Krogstad and Knut-Andreas Lie at SINTEF, as well as Sara Zahedi at KTH for giving me input and helping me during the project.
iii
Abstract
iv
Sammanfattning
Fast marching och fast sweeping för optimala vägval
Contents
1 Introduction 1 1.1 Background . . . 2 1.2 Orienteering . . . 3 1.3 Research Question . . . 5 2 Theory 7 2.1 Formulation of Interface Propagation . . . 72.2 Viscosity Solutions . . . 8
2.3 Hyperbolic Conservation Laws . . . 10
2.3.1 The Linear Wave Equation. . . 10
2.3.2 Higher Order Schemes and Lax-Wendroff. . . 11
2.3.3 The Non Linear Wave Equation . . . 12
2.3.4 Shock-Solution . . . 13
2.4 Basic Algorithm for the Boundary Value Problem . . . 13
2.5 Alternative Mathematical Problem . . . 15
3 Fast Marching Method and Fast Sweeping Method 18 3.1 The Fast Marching Method . . . 18
3.1.1 Casuality Iteration . . . 18
3.1.2 The Update Procedure for the Fast Marching Method . . . 19
3.2 Fast Sweeping Method . . . 20
3.2.1 The Fast Sweeping Algorithm . . . 21
3.2.2 The Iterative Procedure For The Fast Sweeping Method . . . 22
3.2.3 Motivation for the Fast Sweeping Method. . . 23
3.3 Anisotropic Algorithms. . . 25
4 Results 27 4.1 Results from the Fast Marching Method . . . 27
4.1.1 Geodesic Path . . . 29
4.3 Implementation in MRST . . . 31
4.4 Numerical Experiment . . . 33
4.5 Fast Sweeping Route Choices . . . 37
5 Discussion 40
6 Conclusion 43
Chapter 1
Introduction
Fast methods for level sets arise from the studies of interface propagation. The obvious occur-rences of interface propagation include ocean waves, tracing of moving fluid, burning flames, etc. The propagation of an interface can be regarded from an initial value formulation or a boundary value formulation - Level Set Methods are related to the first formulation, and the Fast Marching Method and Fast Sweeping are related to the latter formulation. The problem of propagating interfaces are the foundation of many areas of research, and some applications of level set methods and fast marching methods can be found in separate chapters in [12];
• Level Set Methods in the geometric evolution of curves and surfaces • Level Set Methods applied to the problem of grid generation
• Level Set Methods applied to image enhancement and noise removal
• Level Set Methods/Fast Marching Methods applied to computer Vision regarding shape detection and recognition
• Level Set Methods applied to interface methods for combustion, solidification, fluid me-chanics and electromigration
• Level Set Methods/Fast Marching Methods for computational geometry and computer-aided design
• Fast Marching Method applied to optimality and first arrivals
• Level Set Methods/Fast Marching Methods applied to etching and decomposition in mi-crochip fabrication
In this thesis we will look at optimal route choice for orienteering as an application for the mathematical expression of propagating interfaces. The goal of the thesis is to implement
isotropic and anisotropic algorithms for solving the travel time equation which arise in the re-search of front propagation.
The fast marching method will serve as an algorithm for introducing the isotropic Eikonal travel time equation
F |rT | = 1. (1.1)
This equation will be derived in the next chapter, but as a short summary; T is the arrival time for an expanding front as it crosses each point in a coordinate system and F is the speed function which depends on different arguments related to the physic properties of the problem.
Later in this thesis the focus will be on the fast sweeping method which could be applied to more general anisotropic travel time equations, and to be used in optimal path planning for orienteering.
To the best of author’s knowledge, previously work on the subject of an algorithm for optimal path planning in orienteering are [1], [16] and [3]. In [16] the author addresses the problem of estimating the speed through symbols for the details which an orienteering map consists of, based on gps-data and uses artificial neural networks. This estimation of speed is intended to be used in the calculations of optimal route choices based on an optimal path algorithm. [16] is most of interest regarding further expansion of the work in this thesis. In [1] the author presents a discrete path model approach, and also in [3] a discrete model is used.
The choice of method for this article is numerical experimentation and algorithm construc-tion in Matlab. The use of the MRST toolbox [5] have proven to be a very helpful tool in the implementation in Matlab. OCAD 12 has been used to export symbols from orienteering maps to Matlab, and LiDAR-data has been used in combination with OCAD to produce a digital ele-vation model for the geographical are which the experiments inchapter 4are based on.
The use of general anisotropic algorithms for optimal path planning applied to a complete orienteering map with all symbols would be the subject for future work. In this thesis only some few symbols from an actual orienteering map are used to model the running speed over the computational domain. These symbols are from the author’s point of view the most crucial symbols for deciding a route choice in orienteering. By adding more and more symbols to the model which is presented in this thesis, the model will be more realistic.
1.1 Background
CHAPTER 1. INTRODUCTION 3 functions. Further, these speed functions can be used in algorithms for solving the problem of finding the optimal route choice. If the model used to visualize the mapped terrain is on a highly realistic level, then these algorithms could be an useful tool to do priori analysis on route choices before the gps-data shows which route choice actually is the fastest one. While other well known techniques like Dijkstra’s algorithm or the extended version, A*, are build on a discrete frame-work, the fast methods presented in this thesis are based on a continuous framework.
The basis for the formulation of the optimal path planning problem is the Eikonal Equation. The solution to the Eikonal problem is the first arrival- or viscosity solution. In the context of optimal path planning, other related problems are;
• computation of first arrivals in seismic travel times • construction shortest geodesic paths
• control, computer visibility • line-of-sight evaluations.
The viscosity solution to the Eikonal equation will produce the solution that corresponds to the first arrival of information from the initial disturbance. This means that we can track the solution back to the initial data via the shortest path.
1.2 Orienteering
What follows is a short introduction of orienteering.
orienteering is to draw courses and legs (course between two control points) and analyze what might be the fastest route choices. The fastest route choices are always a combination of the route with the fastest runability and the route choice with the least amount of climb. The com-petitors does not know where the control points are positioned in advance. So all preparations prior to competitions are based on guessing where the controls are positioned and based on the guessed positions, analyze which route choice is the fastest between the controls.
To illustrate the problem each orienteering athlete faces on a typical long leg in an orien-teering course, we look at the 18th leg at the Norwegian National Championship in Ultra Long distance, 20171. The course was placed in the mountains on the west coast of Norway at Kvam-skogen. InFigure 1.1two possible route choices are shown, and inFigure 1.2the elevation pro-file for each route choice are shown.
Figure 1.1: Leg 17 from the Norwegian National Championship in Ultra Long Distance 2017. As we can see fromFigure 1.2the red route choice inFigure 1.1is flat in the beginning, and in the end of the leg there is a big climb. The blue route choice however are shorter, and the climb is taken more evenly over the leg. In general, the shortest route choice is not necessarily the fastest. Running straight under the line connecting the two controls will result in climbing 281m, and running a total length of 2160m. The length of the blue route choice is almost one
1The map used and route choices drawn inFigure 1.1are collected fromhttp://3drerun.worldofo.com/2d/
CHAPTER 1. INTRODUCTION 5
Figure 1.2: Profiles for elevation for route choice 1 (left) and route choice 2 (right). Data collected
from Kartdata©,https://hoydedata.no/LaserInnsyn/.
kilometer shorter than the red route choice and in fact the total climb is less than the red route choice.
Another aspect of choosing a route choice is the steepness of the hills. Lets assume route choice A and route choice B are two route choices with equal length and equal total climb. when speaking of running speed in up- and down hills, the steepness has an impact. If A consists of "smooth" hills, and B of steep hills, then A might be faster even though B will consist of more flat meters. This is because the hills in route choice B might be so steep that it will be difficult to cross them.2
1.3 Research Question
In order to concretize the problems which are addressed in this thesis, the research question is presented below. The overall goal with this thesis is to develop an analysis tool for orienteering which can be used prior to competitions. This analysis tool should will be based on the Eikonal equation, or the travel time equation, and contain a mathematical underlying theory and ap-proach. The research question is the following
• Is it possible to develop an analysis tool based on solving the travel time equation for finding
the optimal route choice on a leg in an orienteering course?
Figure 1.3: Route choice running straight. Data collected from Kartdata©, https://
Chapter 2
Theory
The content ofchapter 2is meant to present the relevant theory for the applications later in the project. In the context of evolving interfaces we will only look into theory relevant for the fast marching method and the fast Sweeping Method. As these methods are evolved from the field of theory of Level Set Methods, the reader will be provided references that covers aspects of Level Set Methods that are of interest.
2.1 Formulation of Interface Propagation
The basic interpretation of level sets is a front T propagating with a speed F . In this section we will have this as a place to start from, and the aim of this section is to visit the basic theory for the algorithms which will be studied later in this chapter.
Consider a boundary that could be a curve T (s, t) in two dimensions or a surface in three di-mensions. The curve T will consist of two regions; a region inside the curve and a region outside the curve, as shown inFigure 2.1. For now, lets consider we have T (0) which is a closed, smooth initial curve. Let T(t) be the family of curves, given by one parameter, which is generated by moving the initial curve T (0) along its vector field with a direction normal to itself with a speed function F (K ) for some properties K of the curve. We want to track the motion as this interface evolves in its normal direction. The speed F (K ) may depend on many factors but for now lets take the influencing factors as local properties L, such as local geometric information, global properties G, which might be properties of the front that depend on the shape and position of the front and independent properties I , which are independent of the shape of the front. Thus lets consider the speed to be a function of L, G and I , i.e. F (L,G, I ).
To track the evolution of the interface, the speed function F must be properly chosen, but at this stage we consider F to be constant and positive. By saying that F > 0 we consider the case that the front moves outwards, or we could say it expands from its initial position. If we consider a two dimensional plane, then we can characterize the position of the front as the arrival time of
Figure 2.1: Propagation of a curve with speed F in normal direction.
the front T (x, y) as it crosses each point (x, y). By using the simple physical relation di st ance =
speed · time we arrive at the equation
1 = FdTd x, (2.1)
Expanding this idea to higher dimensions, the divergence of the arrival times, rT , is orthog-onal to the level sets of T, and we arrive at the boundary value formulation
|rT |F = 1, T = 0 on °, (2.2)
where ° is the initial position of the interface. In the topic of optimal path planning, one may assume that the speed function F depends only on position. Think of planning a route when hiking in the mountains - one plans the route by avoiding obstacles like rivers, lakes, cliffs, etc. One does not consider the local weather in a small valley to have an impact on the route choice. By this assumption that F only depends on position,Equation 2.2reduces to the Eikonal equation. By assuming F also depends on the direction of motion the problem is said to be anisotropic which will be discussed later in this thesis.
2.2 Viscosity Solutions
CHAPTER 2. THEORY 9 We consider the Eikonal equation |rT |F = 1. If the speed F only depends on position x, we are working with a particular case of the general Hamilton-Jacobi equation
Æut+ H(@u, x) = 0 (2.3)
where u 2 R, and H(@u,x) = F |ru| ° (1 ° Æ) (Æ is either zero or one) is the Hamiltonian . We assume that the Hamiltonian H is a smooth function by @u and x. We want to admit non-smooth solutions too to formulate a weak solution to the PDE. One way to do this is to add a viscosity term ≤¢u;
Æut+ H(@u, x) = ≤¢u, ≤ > 0. (2.4)
Now, given a solution u≤we want to show that this is a smooth solution as ≤ goes to zero. Crandall, Evans and Lion [8] defined the weak solution to the Hamilton-Jacobi equation as
Definition 2.2.1. A function u is said to be a viscosity solution to the genereal Hamilton-Jacobi
equation if, for all test functions v:
1. if u ° v has local maximum at a point (x0, t0), then
vt(x0, t0) + H(@v(x0, t0), x0) ∑ 0 (2.5)
2. if u ° v has local minimum at a point (x0, t0), then
vt(x0, t0) + H(@v(x0, t0), x0) ∏ 0. (2.6)
The following practical facts about Definition2.2.1can be found in [8], [7], [2]:
• If u is a smooth solution of the Hamilton-Jacobi equation, then it is a viscosity solution. • If a viscosity solution u is differentiable at some point, then it satisfies the Hamilton-Jacobi
equation at that point.
• The viscosity solution is unique, given appropriate initial conditions.
• The solution produced by taking the limit of the smooth solution u≤as ≤ goes to zero is the viscosity solution.
2.3 Hyperbolic Conservation Laws
As we now have established some knowledge about viscosity solutions, we want to construct correct numerical schemes, and select the correct weak solution corresponding to the viscosity limit of the boundary value equation,
F |rT | = 1.
An extensive resource to this section can be found in LeVeque [4].
The study of hyperbolic conservation laws are important to be able to pick the correct nu-merical solution of hyperbolic differential equations. The characteristics of the hyperbolic dif-ferential equation will play an important role in this part, and it is necessary to know how the characteristics influence the choice of numerical schemes for the problem.
2.3.1 The Linear Wave Equation
An example of a hyperbolic conservation law, is the linear one dimensional wave equation
ut+ ux= 0, u(x,0) = f (x), t > 0 and x 2 R. (2.7)
The exact solution toEquation 2.7is u(x, t) = f (x ° t) i.e. a convection wave that starts from its initial position and convects with the same form as t increases. We now want to develop a numerical scheme for this equation. We start by dividing the x-t space into a grid with space step ¢x and time step ¢t. Lets denote the numerical solution for index i at time n as uni.Depending on which direction we discretize using Taylor expansion, we arrive for the time derivative at the schemes, • Forward difference: ut=u n+1 i °uni ¢t +O(¢t) • Backward difference: ut =u n i°un°1i ¢t +O(¢t)
• Central difference: ut=un+1i °un°1i
2¢t +O(¢t).
The same discretization can be done in space to obtain • Forward difference: ux=ui +1n °uin
¢x +O(¢x)
• Backward difference: ux=uni°ui °1n
¢x +O(¢x)
• Central difference: ux=uni +1°uni °1
CHAPTER 2. THEORY 11
Figure 2.2: Solution toEquation 2.7, which is constant along the slopes.
The solution toEquation 2.7, u(x, t) = f (x ° t), means that the solution u at any point x at time t is the same as the value of the initial data at the point x ° t on the x-axis. To illustrate the solution, we can say that the solution u is constant on lines with slope 1 in the x-t plane. So as time increases, the solution at point A inFigure 2.2is the same as the solution in point B. We say that the domain of dependence of point B is the point A.
By the principle of domain of dependence, the forward difference will not capture the solu-tion in a correct manner. But the backward difference is preferable since it uses values upwind of the information propagation. These straight lines are known as characteristics.
2.3.2 Higher Order Schemes and Lax-Wendroff
By using Taylor expansion it is possible to produce higher order schemes. Taylor expanding
ut = °ux yields, u(x, t + k) = u + kut+k 2 2 ut t+O(k 3) = u ° kux+k2 2 uxx+ o(k 3), (2.8)
and this equation can be formulated (by following the notation in [12]) as the following scheme uin+1= uni ° kD0xuin+k 2 2 D +x°xun i, (2.9) where D0xun
i is shorthand foru(x+¢x,t)°u(x°¢x,t)2¢x and D+x°xuinis shorthand for u(x+¢x,t)°2u(x,t)+u(x°¢x,t)
2.3.3 The Non Linear Wave Equation
As the characteristic solution is an important concept of the fast marching method and the fast sweeping method, we will in this section look at a "model problem" for illustrating the principles of how to use the characteristics. This model problem is the non-linear wave equation.
If the wave speed a depends on the position x, i.e.,
ut+ a(x)ux= 0, u(x,0) = f (x), a(x) 6= 0, (2.10)
we have the non linear wave equation. When dealing with this equation, the choice of nu-merical schemes should depend on the sign of a. If a > 0 the information travels from left to right, and a backward difference operator is preferred, i.e.
D°xu ¥u
n i ° uni °1
¢x . (2.11)
If the wave speed is negative, i.e. a ∑ 0, the information travels from right to left, hence the correct difference operator to use is the forward difference operator,
D+xu ¥u
n i +1° uni
CHAPTER 2. THEORY 13
2.3.4 Shock-Solution
By saying that the wave speed depends on the solution u itself, we are studying the non-linear "Burger’s Equation", and in this case we look at the "Inviscid Burger’s Equation",
ut+ uux= 0, u(x,0) = f (x). (2.13)
If we consider a particle moving in the x-t plane, and the movement is parameterized by s, we get by using the chain rule,
du(x(s), t(s)) d s = ux d x d s + ut d t d s = d t d s ∑ ut+d x d tux ∏ . (2.14)
If we suppose that the trajectory of the particle is such thatd xd t = u we identifyEquation 2.13
thus the right hand side is zero and dud s = 0. Thus u is constant along the trajectories which is known as "characteristics". Furthermore the slope of these trajectories does not change, hence these characteristics are straight lines.
However, the characteristics does not need to be parallel. For example in one part of the domain, lets say x = [0 1] the characteristics are lines with slope one, i.e. 45 degrees, and in
x = [1 2] the characteristics are vertical lines. At some point in time, the slopes will cross each
other, and when this happens, shocks will appear. If the characteristics are spreading out - for example for x = [0 1] the characteristics are pointing 45 degrees to the left, and for x = [1 2] the characteristics are pointing 45 degrees to the right. The solution in this region will be charac-teristics that are spread out like a fane. This type of solution is referred to as a rarefaction wave. Illustrations and more detailed explanation can be found in [4].
2.4 Basic Algorithm for the Boundary Value Problem
Up to this part of this chapter we have focused on different theoretical principals which the solution to the travel time equation depends on. In this section the numerical scheme which solves the travel time equation, and scheme which the fast marching method is based on will be presented.
We start by recalling the boundary value problem
F |rT | = 1. (2.15)
By introducing the general Hamilton-Jacobi equation
ÆUt+ H(Ux,Uy,Uz, x, y, z) = 0 (2.16)
H(Ux,Uy,Uz, x, y, z) = F q
Ux2+Uy2+Uz2° 1. (2.17)
If we focus on the general Hamilton-Jacobi equation in one dimension we have the form,
ÆUt+ H(Ux) = 0 (2.18)
If we set Ux= u and differentiate we get the hyperbolic conservation law,
Ut+ [H(u)]x= 0. (2.19)
H(u), or in the general cases G(u), is defined as the flux function, and it is clearly shown in
[4] and [12] that numerical fluxes g can be constructed to approximate the above equation by,
ui +1n ° uni
¢t =
g (uin,ui +1n ) ° g(ui °1n ,uni)
¢x . (2.20)
As described insubsection 2.3.3the choice of difference operator should depend on which
direction the information travels. We saw that for positive wave speed a > 0, the backward dif-ference D°xuni was chosen. When the direction of the wave speed is not entirely known at all times, an adaptive construction of a difference scheme can be chosen. This adaptive construc-tion of a difference scheme will ensure that the scheme only evaluates informaconstruc-tion propagating upstream. Recall the non-linear wave equation with wave speed a(x),
ut+ a(x)ux= 0, u(x,0) = f (x) and a(x) 6= 0.
Now we can choose the correct direction of the upwinding so that the difference scheme always includes the mathematical domain of dependence, by writing the update scheme,
uin+1= uni ° ¢t[max(0, ai)D°xuni + min(0, ai)D+xuin]. (2.21) We see from the update scheme that when a > 0 the expression max(0,ai) = aiand min(0, ai) = 0 thus the difference operator D°xun
i is used. And, consequently, when a < 0 the expression max(0, ai) = 0 and min(0,ai) = ai thus the difference operator D+xun
i is used. This idea can be used to rewrite theEquation 2.2as,
2 6 6 4 max(D°xun i j k,0)2+ min(D+xuni j k,0)2 max(D°yun i j k,0)2+ min(D+yuni j k,0)2 max(D°zun i j k,0)2+ min(D+zuni j k,0)2 3 7 7 5 1 2 =F1 i j k (2.22)
CHAPTER 2. THEORY 15
2.5 Alternative Mathematical Problem
Before continuing with the theory specific to the fast marching method and the fast sweep-ing method we will first look at an alternative description of the travel time equation we study. Until now we have spoken about isotropic speed functions related to the travel time equation. This means that the cost of moving in any arbitrary direction costs the same. In this section anisotropy related to movement in the domain will be introduced.
As an example of anisotropy we can consider a road in a steep hill side. Often when roads are constructed in steep hill sides, the road are constructed as serpentines. This is because it is easier to walk upwards in serpentines instead of straight up the hill side. But if we consider the case when we want to walk down the hill side, we often don’t care about the road, and just sets off straight down the hill side. The same principle is often used in theory related to some kind of navigation through a cost matrix - and if the speed depends of which direction we move we call the speed function anisotropic.
In general the fast marching method and the fast sweeping method will not be applicable to anisotropic metrics. The details will be described in the next chapter where the specific theory are described in detail.
Definition 2.5.1. When the cost W (i , j ) of moving from one grid point xi ,jto another grid point
xs,t depends on the current position, namely xi ,j applied to some reference system, the metric
W is said to be isotropic.
Definition 2.5.2. When the cost W (i , j ) of moving from one grid point xi ,jto another grid point
xs,t depends on the current position and the direction of movement, the metric W is said to be anisotropic.
As an alternative for the travel time equation mentioned previously in this thesis, we will now derive a slightly different travel time equation.
Consider a given surface z = S( ¯x) as the one inFigure 2.3at a grid point ¯x0. In this grid point
the steepest direction of ascent is defined as rs
|rs|= ¯q (2.23)
Assume we have the following velocities moving upwards-, downwards- and transversal di-rection,
• vu, upward direction along q,
Figure 2.3: Example of contours on an orienteering map.
If we now restrict ourselves to the neighbourhood around the grid point ¯x0with constant rs,
a possible travel-time function has the form,
T ( ¯x) =q( ¯x ° ¯x0)TK°2( ¯x ° ¯x0) ° g( ¯x ° ¯x0)Tq,¯ (2.24) where, • K = Q √ a 0 0 b ! QT, Q = [ ¯q, ¯q?], • a =≥12≥v1u+v1d ¥¥°1 , b = vt and • g =12≥v1d°v1u ¥ (gravitation term)
We now have the framework for the alternative travel time equation. Taking the gradient of
T gives, rT = p 1 ( ¯x ° ¯x0)TK°2( ¯x ° ¯x0) K°2( ¯x ° ¯x0) ° g ¯q = 1 |K°1( ¯x ° ¯x0)||° g ¯q (2.25)
By moving the gravitation term to the left side of the equation it follows that,
CHAPTER 2. THEORY 17 Now the velocities vu, vd and vt are given in each grid point ¯x.
Remark: By defining the speed function F as F = vu = vd = vt the problem reduces to the isotropic Eikonal equation,
F |rT | = 1. (2.27)
Fast Marching Method and Fast Sweeping
Method
3.1 The Fast Marching Method
Before we start off this section, a maybe confusing issue is visited. The Eikonal equation is a boundary value problem. However the scheme we are going to build for the fast marching method could be interpreted as an initial value problem. The reason for this is that we are go-ing to construct the method where the solution is built outward from the boundary value data hence the initial value interpretation. Anyway this is only the algorithmic construction that has this propagation view.
The study of fast marching method will be continue to be studied in this chapter and we will also introduce the fast sweeping method. The fast sweeping method will be the most important algorithm of the two, it will be used in chapterchapter 4to generate the optimal route choices on an orienteering map.
3.1.1 Casuality Iteration
This chapter will start with an important concept for the class of single-pass algorithms, which the fast marching method is, namely causality. Recall equation (2.22). By following the upwind scheme given in [11] we rewrite equation (2.22) as
2 6 6 4 max(D°xTn i j k,°D+xTi j kn ,0)2 +max(D°yTi j kn ,°D+yTn i j k,0)2 +max(D°zTi j kn ,°D+zTi j kn ,0)2 3 7 7 5 1 2 =F1 i j k. (3.1)
where still T is the travel time for a propagating front.
CHAPTER 3. FAST MARCHING METHOD AND FAST SWEEPING METHOD 19 To solve equation (3.1), an approach is to consider all the neighbor grid points around Ti ,j,k, i.e. Ti +1,j,k, Ti °1,j,k, Ti ,j +1,k, Ti ,j °1,k, Ti ,j,k+1, Ti ,j,k°1. We assume all the neighbor points are known, so we solve equation (3.1) for Ti ,j,k. Assuming we have N2grid points in the computa-tional domain, and the fact that we need some iterations to obtain a solution with a satisfactory approximation, we can construct the loop as,
for i t = 1 :n
for i , j , k = 1 :n
compute T from the Quadratic equation given a l l neighbours end
end
where the quadratic equation is equation3.1.
This approach does not take advantage of the central idea behind the Fast Marching Method which is to use upwind values to construct the solution Ti ,j,k. The essence behind Fast March-ing Method is to solve equation (3.1) from the smallest T value, and move from the boundary solution using upwind values. This means that we don’t need to iterate to obtain a satisfactory approximation. The way we move from the boundary value and towards the starting position is explained here; we start at the smallest T and then,
1. update solution downwind 2. compute new possible values
3. choose the nearest grid point with the smallest value 4. freeze this point and continue with 1.
As we move out from the boundary value grid point, the upwind structure ensures that we will not move to a grid point containing a larger value of T than the previous evaluated grid points. We are marching the solution outwards and visits all the points with minimum trial value for T . This way of computing the solution is obviously a very large difference to Dijkstra’s graph algorithm which revisits already evaluated nodes to see if there is a cheaper edge to move along.
3.1.2 The Update Procedure for the Fast Marching Method
start_point = S T r i a l [ ] = a l l _ g r i d _ p o i n t s while A != start_point A = min( T r i a l [ ] ) T r i a l [ ] . pop(A) Known [ ] . append(A)
for a l l neighbors N to A not in Known[ ] i f N in Far [ ] T r i a l [ ] . append(N) else computeNewT(N) end end end
The restriction for how fast this method works, is how to obtain the grid points with the lowest T value. An example is to use the min-heap structure presented in [12].
3.2 Fast Sweeping Method
The fast sweeping method, which is the mainly applied algorithm in this thesis, will be described in this section. The fast sweeping method is more interesting to use in the evaluation of optimal route choices in orienteering because it has the ability to handle anisotropy on a more gener-ally basis than the fast marching method. As the construction of a proper speed function from an orienteering map will consist of high anisotropies, the fast sweeping method will be better suited for the problem which is studied in this thesis.
Recall the causality approach in the fast marching method where the update of the solution follows the causality in a sequential way. As described above the solution is updated one grid point at a time in the order that the solution is strictly increasing1.
The Fast Sweeping Method is an iterative algorithm and was first used to compute the
dis-tance function in [15]. The Fast Sweeping Method uses nonlinear upwind differences in
com-bination of Gauss-Seidel iterations with alternating sweeping ordering. Compared to the fast marching method, the fast sweeping method follows the causality along characteristics in a par-allel way. The characteristics are grouped according to their directions, where each Gauss-Seidel iteration with a specific sweeping ordering covers a group of characteristics simultaneously. For
N2discretization points in the computational domain, the Fast Sweeping Method has a
compu-tational complexity of O(N) at it’s best.
CHAPTER 3. FAST MARCHING METHOD AND FAST SWEEPING METHOD 21
3.2.1 The Fast Sweeping Algorithm
The following sections follows the theory presented in [15]. We consider the model problem,
|ru|F = f , T = 0 on ° (3.2)
where f (x) > 0. xi ,j is used to denote a grid point in the computational domain ≠, h as grid size and ui ,jh denotes the numerical solution at xi ,j. To discretize the partial differential equation at interior grid points, we use the Godunov upwind scheme proposed in [14] and the discretization becomes for the interior points i = 2,..., I ° 1, j = 2,..., J ° 1,
[(uhi ,j° uhxmi n)+]2+ [(uhi ,j° uhymi n)+]2= fi ,j2 hh (3.3) where in the x-direction uhxmi n= min(ui °1,jh ,uhi +1,j), in the y-direction uhymi n= min(uhi ,j °1,ui ,j +1h ) and (x)+=
8 < :
x, if x > 0
0, if x < 0. At the boundary of the domain we use an one sided difference.
In the same way as in the Fast Marching Method, the solution is initialized by assigning exact values at grid points in ° to enforce u(x) = 0 for x 2 °. At all the other grid points the solution is initialized with large values2.
The sweeping which assigns the correct values at each grid point xi ,j works as follows: at each xi ,j who’s value is not fixed during the initialization we compute ¯u of euqation (3.3) using the values of uhi ±1,j and uhi ,j ±1followed by updating ui ,jh according to uupd atedi ,j = min(ui ,jold, ¯u).
In two dimensions we sweep this update routine according to the four orderings 1. i = 1 : I, j = 1 : I
2. i = I : 1, j = 1 : I 3. i = I : 1, j = I : 1 4. i = 1 : I, j = I : 1
The unique solution to equation3.3is
¯ u = 8 > < > :
mi n(uhxmi n,uhymi n) + fi ,jh |uhxmi n° uhymi n| ∏ fi ,jh uh
xmi n+uhymi n+
q 2f2
i ,jh2°(uhxmi n°uhymi n)2
2 |uhxmi n° uhymi n| < fi ,jh
. (3.4)
For the solution in n dimensions we refer to the procedure presented in [15].
3.2.2 The Iterative Procedure For The Fast Sweeping Method
We consider the same computational domain as for the Fast Marching Method. In two dimen-sions the direction of sweepings could be as stated above,
1. i = 1 : I, j = 1 : I 2. i = I : 1, j = 1 : I 3. i = I : 1, j = I : 1 4. i = 1 : I, j = I : 1
and the algorithm for the Fast Sweeping Method would be, sweep_xdir = [ 1 :Nx ; 1 : Nx ; Nx: °1:1;Nx: °1:1];
sweep_ydir = [ 1 :Ny;Ny: °1:1;1:Ny;Ny: °1:1]; for each i in sweep_xdir
for each j in sweep_ydir
i f u( i , j ) not on boudary
uxmin = min(u( i °1, j ) ,u( i +1 , j ) ) uymin = min(u( i , j °1) ,u( i , j +1)) else
uxmin = u( c l o s e s t to x°boundary ) uymin = u( c l o s e s t to y°boundary ) end
i f | uxmin°uymin| >=W( i , j )*h
ubar = min(uxmin , uymin)+W( i , j ) *h else
ubar = (uxmin+uymin+sqrt (2*W( i , j )^2*h^2°(uxmin°uymin)^2))/2) end
u( i , j ) = min(u( i , j ) , ubar ) end
end
CHAPTER 3. FAST MARCHING METHOD AND FAST SWEEPING METHOD 23
3.2.3 Motivation for the Fast Sweeping Method
As previously mentioned the Godunov upwind difference scheme enforces the causality. The boundary conditions enforces the propagation of information to move outwards from ° since
°is contained in the computational domain. The idea behind the Gauss-Seidel iterations with
different alternating sweeping orders, is that each sweep follows the causality of groups of char-acteristics in certain directions simultaneously. The update rule will only update the solution at grid points which has not reached it’s minimum value thus all other grid points will not be changed. When a grid point finally reaches it’s minimum value, the value will not change in later iterations.
To illustrate the motivation of the fast sweeping method we consider the eikonal equation in one dimension,
|rd(x)| = 1, d(x) = 0, x 2 °. (3.5)
In this example the characteristics follows straight lines which radiates from the set °. Since this problem is one-dimensional only two sweeps are needed to obtain the solution, i.e. from left to right, and then sweep from right to left. By sweeping in the two directions we have captured the solution at the two possible directions of the characteristics. More accurately
speaking, the distance function at a given point xi can be computed from it’s two neighbours
by di = min(di °1,di +1) + h. The first sweep will cover the right-moving characteristics, and the second will cover the left-moving characteristics.
In higher dimensions the case is more delicate as the characteristics may have infinitely many directions, which may not be captured exactly by a cartesian grid. In this case we need to address how many Gauss-Seidel iterations are needed to obtain a solution, and how large is the error for this solution.
Figure 3.1: Figure showing the four different groups of characteristics; up-right, up-left, down-left and down-right.
In this example the right hand side cost function f is set to one in all coordinates, which means that there is no extra cost involved to move through one grid cell xi ,j compared to an-other. This means that the characteristics are “simple”. For more general cost function one may need more than one sweep to cover one characteristic curve. By following the motivation above, and by considering the equation above, we can predict the fast sweeping method to converge
after 2n sweeps for problems in n-dimensions. Now, as we are speaking of more general cost
functions, we might need more than 2nsweepings to cover the characteristics.
The order of dependence of the grid points are shown in figure3.2. The red dots represent interior grid points of the domain while the blue dots represents grid points which might be influenced by boundary conditions. The black line marked with “1” shows that the two grid points connected with the line are dependent of the left-down grid point. In the same way we see how the grid points connected with the black line marked with “2” are dependent of the two previously computed grid point values. This figure serve as a general illustration of how the solution is computed over the domain as the algorithm is swept over the grid points by varying
CHAPTER 3. FAST MARCHING METHOD AND FAST SWEEPING METHOD 25
Figure 3.2: Figure showing the order of dependence of the computational domain. The figure specifically shows the order of dependence for the upper-right part of the domain.
3.3 Anisotropic Algorithms
Recall that if the cost of moving away from a grid point only depends on the grid point, and not the direction of motion, the travel time equation we wish to solve is isotropic. Hence the optimal path will be given by solving the travel time equation and move in the direction of steepest descent of the travel time T . The fast marching method and the fast sweeping sethod has been the focus in this thesis. While the general fast marching method is relying on one direction of the characteristics, the fast sweeping method can handle multiple directions of the characteristics by sweeping over the domain until all the directions of characteristics are treated.
In the update procedure of the fast marching method described previously in this chapter,
subsection 3.1.1, the updated value of a grid point must be based on information of neighbour-ing grid points. These neighbourneighbour-ing grid points must be of smaller value in order to compute
¯
ui ,j in increasing order. The violation of this causality property is the main reason for why the fast marching method will fail in the general case of an anisotropic travel time equation. In the case of an anisotropic travel time equation, the characteristics will not be straight lines radiating from the boundary.
The Ordered Upwind Method [13] is a promising start for an algorithm to solve a general anisotropic travel time equation applied on an orienteering map. This method relies on a trian-gulated mesh. Compared to the fast marching method and the fast sweeping method, the or-dered upwind method is much more intricate to implement. Where the fast marching method could violate the causality principle, the ordered upwind method searches for smaller values from neighbouring grid points along an active front. This search along the active front has the purpose to find a set of neighbouring grid points (which not necessarily needs to be a direct neighbour of the current grid point xi ,j) whose value have been accepted, and then construct a simplex with these grid points to update the numerical solution T . The Ordered Upwind Method computes the numerical solution T to the general anisotropic travel time equation by evaluating subsets of the computational domain. The Ordered Upwind Method includes the three subsets;
• a set where T is known,
• a set where tentative values for T has been computed, • a set where no information regarding T is known.
Chapter 4
Results
In this chapter the results will be presented. First a more detailed description of the implemen-tation will be described, and then some results from the compuimplemen-tation will be presented. The first part of this chapter will present some results from the fast marching method. As the fast marching method is the main motivation for the ordered upwind method, some results from the implementation of the fast marching method is presented in this first part of the chapter.
4.1 Results from the Fast Marching Method
The goal of the implementation of the Fast Marching Method and Fast Sweeping Method was to, • have an interesting example to execute the algorithm on,
• being able to solve the travel time equation F |rT | = W and
• solve the travel time equation with isotropic- and anisotropic speed functions.
The test-data which the fast marching method in this part considers are images. These pictures are rgb-colored pictures which where converted to a grey-scale-valued matrix which worked as the metr i cs or cost ° matr i x for crossing a pixel, i.e. W (x) f or x 2 ≠ where ≠ could be any two-dimensional sized domain. The aim for the code was to have a designated start point x0, and an end point x1, and then the methods would march/sweep over the domain
to obtain the shortest path ∞,
min ∞ Z ≠W (∞(¬)| d∞ d¬|d¬. (4.1)
The values in the cost-matrix W (x) depends on where we define x1. It takes the value from
the pixel which x0represents, and then all the other values in W (x) are given by taking the
dif-ference between the pixel value in x0and the pixel value in all the other points in ≠.
Defining the result from the Fast Marching as the "shortest path" might be counter-intuitive since the shortest way from one point to another is always a straight line between these points. But in this context "shortest path" is referred to as the the path ∞ which minimizes the weighted length in relation to the cost-matrix.
(a) Initial image. (b) Grey scale image.
(c) Image with starting- and end point.
Figure 4.1: (a):image showing a system of capillary blood vessels. The dark areas are tissue which is hard for the blood to move through and light areas are easy to move through. (b): The same vessels asFigure 4.1abut as grey scale image.(c): Start- and end point defined. The initial image
is downloaded fromhttps://goo.gl/cUgWfr.
CHAPTER 4. RESULTS 29
4.1.1 Geodesic Path
As described inchapter 2, the fast marching method will start from the arrival point1and
cal-culate the arrival time T from the travel time equation for every point in the domain. The result from solving for T can be displayed as a distance map and a geodesic map.Figure 4.2shows the result from computing the distance map and the gradient map.
(a) Distance map. (b) Gradient map of the distance map.
Figure 4.2
The distance map,Figure 4.2a, has different colors which represents the arrival times to x1.
Blue means short arrival time and dark red means long arrival time. So the fastest way is to move as much as possible through fast colors and avoid red areas.
To extract the geodesic curve we use a numerical approach to perform a gradient decent. We compute the geodesic curve between any points in the domain through the gradient decent,
∞0(¬) = °ø¬rD(∞(¬)) (4.2)
where ∞(0) = x1and ø¬> 0 is the step size we use to iterate the curve. rD(∞(¬) is the gradient
ofFigure 4.2a. We obtain a discrete curve,
∞k+1= ∞k° øG, (4.3)
where G = rD(∞(¬) and ∞kis the approximation ofEquation 4.1.1at ¬ = øk.
FromFigure 4.3we see the geodesic curve that connects the departure- and arrival coor-dinate, respectively x0and x1. The geodesic curve does not move perfectly through the large,
Figure 4.3: Figure showing the geodesic path between x0and x1.
white vessel. By comparing Figure 4.2a and Figure 4.2b to Figure 4.3, we see that the
com-puted geodesic curve follows the areas in the distance map which has most blue coloring; by approaching the red dot from right, we cross less yellow- and green colored pixels compared to an approach from left, and by considering the area around the green dot; we see that taking a right route choice involves crossing more dark blue colored pixels than taking a left route choice.
4.2 A Comparison of the Fast Sweeping Method and Fast
March-ing Method
We have now had a look at the results from the fast marching method and in this section we do a comparison with the fast sweeping method with the right hand side metric given by,
W = 1 + 8e°x2+y22æ2 . (4.4)
CHAPTER 4. RESULTS 31 50 100 150 200 250 300 50 100 150 200 250 300 FMM
(a) Generated route choice by using the Fast March-ing Method. 50 100 150 200 250 300 50 100 150 200 250 300 FSM
(b) Generated route choice by using the Fast Sweep-ing method.
The rest of the examples in this chapter will be based on anisotropic metrics, thus the fast marching method will not be studied any further in this chapter. The rest of the examples are solutions based on the fast sweeping method, and will show the results of the geodesic paths between two grid points on an orienteering map.
4.3 Implementation in MRST
In MRST the Finite Volume Method is used to solve the system of equations. MRST is a free open-source software for reservoir modelling and simulation, hence some terms related to the physics of reservoir theory is used; "rock structure" which is a representation of a porous medium and; "permeability" which is defined as a porous medium’s ability to transmit fluids. Regarding this thesis and the implementation below, the permeability is analogous to the speed of moving in a certain direction at a grid point xi j. In the implementation using MRST the following key tools are used,
• To set up the computational domain we define a cartesian grid defined by the number of discretization points N x and N y, as well as the scaling of the cells, namely Dx and
D y. The geometry is generated as a struct consisting of cells, faces, nodes, and general
specifications about the grid.
the grid points on the surface is then calculated in x- and y-direction. After this step, we use the mean values of the "height" on the surface to compute the gradient of the surface rS.
• After the surface is defined we move over to defining the speed function F . The default value is defined to vf l at = 3.5m/s, which will work as a basis for the speed. This is the speed which is adjusted accordingly for movement up-, down- and transversal to con-tours. To account for moving in the upward, downward and transversal direction we de-fine different pace-multiplicators which adjust vf l at. InTable 4.1the pace-multiplicators are tabulated. Note that pt are defined as the average of puand pd.
G0 0.0 0.1 0.15 1
pu 1 1.25 1.45 8
pd 1 0.9 0.95 3
pt 1 1.075 1.2 5.5
Table 4.1: Pace multiplicators for the speed according to direction of motion and gradient of surface.
Now as the speed is specified for the tabulated values, the values for the speed in all the grid points in the domain is defined by dividing vf l at with the interpolated value of
vu/vd/vt according to the surface at this grid point.
• When the speed at all grid points are adjusted according to the surface S we define a rock structure to represent the speed function F as permeability in every grid point. From the documentation in [5], the rock property is characterized by the permeability that gives the rock’s (grid cell’s) ability to transmit fluid (in our case the speed for moving from one grid point to a neighbour grid point).
• The next step is to compute the solution. As this method is based on the MRST toolbox, the solution is produced by iterations of an initial solution until the solution has converged.
As a start guess we use the method digraph [9] which takes the computational domain
as an argument for constructing the graph, and then the transmissibilities as values for crossing a link between two nodes in the graph. Then the equation is solved by iterating the solution until a given error tolerance between the iterations is reached.
CHAPTER 4. RESULTS 33 InFigure 4.5we see the contours generated from LiDAR data2, ©Kartverket. The blue con-tours are set as zero-level and the yellow concon-tours have the highest altitude. From the figure we see that the south-east part consists of large white areas which means this area is flat - and in reality some parts of this area is a lake, and rivers which are in connection to the lake.
50 100 150 200 250 300 350 400 50 100 150 200 250 300 350 400
Figure 4.5: Contour plot of the test data obtained from ©Kartverket.
In the areas where there are lakes and rivers which are not passable, we set the running speed to zero, because we want to avoid crossing these areas. InFigure 4.6b,Figure 4.6aandFigure 4.7
we see the velocity related to moving in upward-, downward- and transversal direction. The dark blue areas are impassable. By comparingFigure 4.5andFigure 4.7we see that the area, which corresponds to the impassable areas of water, does not overlap with the zero velocity shown as dark blue inFigure 4.7. This offset can also be found by looking at the position of the river in the valley in the west part ofFigure 4.6a. The reason for the offset is that the symbols from the orienteering map is not geo-referenced exactly3according to the LiDAR data.
In the first two figures we show the travel time from the starting point to the defined end point. The color bar shows the levels of travel time where blue is the color from the start point and at the end point the color is yellow. So by following the contours we see the result which is generated in the figure with the region of fastest route choice.
4.4 Numerical Experiment
So far in this chapter we have been looking at the practical results of the fast sweeping method and the algorithm has been applied to an orienteering map with and without
obstacles/impass-2LiDAR data in Norway is distributed throughhttps://hoydedata.no/LaserInnsyn/
(a) Plot of the downward velocity. Dark blue regions are "impassable regions".
(b) Plot of the upward velocity. Dark blue regions are "impassable regions".
Figure 4.6
able areas. This section is dedicated to look at the numerical performance of the fast sweeping method. While the fast marching method is a single pass algorithm, the fast sweeping method is an iterative algorithm, so to check the convergence while iterating the algorithm is a natural choice of a numerical analysis.
To check the convergence and the numerical performance of the fast sweeping method, the same exact solution as in [15] where used, i.e.,
u(x, y) = Ø Ø Ø Øe ≥p x2+y2°r0¥(ax2+2bx y+c y2+d ) °1ØØØØ (4.5)
where r0= 0.2, a = 0.1, b = 0.5, c = 1 and d = 0.4. Recall the equation we wish to solve,
|ru| = f .
So we use f (x, y) = ru(x, y). The calculations of ru is carried out in the following way,
CHAPTER 4. RESULTS 35
Figure 4.7: Plot of the transversal velocity. Dark blue regions are "impassable regions".
n Error h 50 0.1099 0.04 100 0.0535 0.02 200 0.0264 0.01 400 0.0131 0.005 800 0.0065 0.0025 1600 0.0033 0.00125
Table 4.2: Table or the error which shows first order convergence.
@u @x = √ x(0.1x2+ x y + y2+ 0.4) p x2+ y2 + (x + 2y) µq x2+ y2° 0.2 ∂! · ... ...e ≥p x2+y2°0.2¥(0.1x2+x y+y2+0.4 ) ·... ... µ e ≥p x2+y2°0.2¥(0.1x2+x y+y2+0.4 ) °1∂±µe ≥p x2+y2°0.2¥(0.1x2+x y+y2+0.4 ) °1∂ (4.7)
By taking the square root ofEquation 4.6andEquation 4.7squared we get, ru = s @u @x 2 +@u@y 2 . (4.8)
By usingEquation 4.8 as f (x, y) in the fast sweeping algorithm and by denoting ˆu as the
numerical solution, we compare the numerical solution ˆu to the exact solution u. The error
is calculated as e1= k ˆu°uexactk1
kuexactk1 . Table 4.2 shows the error, and inFigure 4.9a loglog-plot of
(a) Plot of the time from the starting point with co-ordinates (50,0) to the end point with coco-ordinates (400,350). Simulation performed in MRST.
(b) Plot of the generated route. In this simulation the route is visualized through the blue region, so the optimal route is contained in the blue region. Simu-lation performed in MRST.
CHAPTER 4. RESULTS 37
4.5 Fast Sweeping Route Choices
The results in this section are comparable with the results insection 4.3containing the results from the MRST based implementation. The areas of the map with cliffs, lakes, rivers and pri-vate property are still defined as impassable and the velocity of moving through such regions is therefore set to a slow pace. In the results from the computation of the geodesic paths it should be expected to see the path avoid such areas.
As a first look at the result we can see fromFigure 4.10bthe start grid cell where the level fronts starts to evolve from the boundary, and we can see the geodesic path generated from
following the steepest gradient descent. In Figure 4.10a and Figure 4.11a we see two route
choices generated by using the fast sweeping method, whereFigure 4.11ais the one related to
Figure 4.10b. In these figures we see the contours from the terrain model. In these two examples we see two typical scenarios one can meet. InFigure 4.10athe fastest route choice is a straight choice where climbing seems to be faster then following the height and running transversal and follow the contours. InFigure 4.11ait is possible to exploit flat ground in the south east area, and the generated route choice picks a route covering this flat area.
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 FSM Route Choice
(a) Generated route choice using the Fast Sweeping Method for grid size 100x100.
FSM Level Plot 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 FSM Route Choice
(a) Generated route choice using the Fast Sweeping Method for grid size 100x100. In this simulation the Fast Sweeping Method generates a route choice that takes as litle climb as possible.
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100
(b) Level curves from the route choice using the Fast Sweeping Method for grid size 100x100.
In the next example we use the metrics given by,
W = 55max(°0.2,sin(40x ) + cos(y ° x 2 50 )) + 0.02x + y ° 500 50 2 . (4.9)
By adjusting the pace-multiplicators we should expect to see the route choices changes. For
example if we adjust pu in such a way that the velocity of moving upwards is decreased we
should expect the route choice to increasingly follow the hill sides rather than crossing contours and run more straight. By changing the upward velocity pufrom,
• old pu= [1, 1.25, 1.45, 8] • new pu= [5, 3, 5, 10],
CHAPTER 4. RESULTS 39 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 FSM Route Choice
(a) Generated route choice using the Fast Sweeping Method for grid size 100x100, where the old pu is
used. 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 FSM Route Choice
(b) Generated route choice using the Fast Sweeping Method for grid size 100x100, where the new pu is
Discussion
We have seen through the previous chapters that the fast sweeping method can be applied to the problem of picking a route choice in orienteering. The data set which have consisted of a digital elevation model generated from LiDAR data and symbols from the file of an existing orienteering map, have proven to be able to be used as input in the travel time function we have studied. In this chapter we will address some of the issues with the model and what can be done to make the model more realistic.
The question of which route choice is the fastest will always have a complex answer, or maybe it might be more correct to call it a suggestion. Between athletes there are a difference in physical capacity and orienteering technical capacity. Within these two types of athletes there are differences too. It might be natural to assume the athletes in the start field to be equally fit and the differences between the athletes are related to different specialties. Some athletes have a specialty of running in tough terrain areas with slow running speed as marches or uphills, while other athletes have a specialty of running fast on hard surfaces with high running speed such as a path or a road. When a course maker is designing a leg he often wants to challenge the athletes to make a decision between the route choices which utilizes their own specialties. When speaking of orienteering technical abilities some athletes might attack the controls along a technical demanding route choice, or some might take a less technical demanding.
A statistical research of picking the best route choice can be found in [10]. In this article he addresses the question of what really is the best route choice, and a number of different bias is mentioned regarding the correctness of the results from the study. One major bias when studying route choices is the fact that most runners will pick the route choice which suits their personal specialties the best. As Myrvold states: "In some ways the runner’s task is easier than the coach’ when it comes to route choice analysis. The runner needs only to find what will work best for him or her, while the coach is usually looking for some "best" route choice". A perfect model which accounts for these different types of specialties between the runners might be too complex to implement, but some adjustments should be possible to apply on for instance the
CHAPTER 5. DISCUSSION 41 running speed setup in the implementation of the travel time algorithm, such as the values for
vu, vd and vt. In this thesis we have not compared actual running times to the running time generated by the fast sweeping method. A perfect model will hopefully be able to calculate a "superman" time on a given leg given that the runner we want to model picking the optimal route choice is fully rested before running the leg, and is capable of running in full speed with-out slowing down because of exhaustion. However this perfect model might meet many obsta-cles in the implementation. The model needs to account for the difference in running speed crossing the same map symbol, which is a fact as some terrains are more wet than other terrains and some terrains consist of more under-vegetation which might not be accounted for in the mapping e.g.
As mentioned throughout this thesis, the data which the speed function F depends on is based on a digital elevation model and some impassable symbols. Important symbols which have not been a part of the speed function F are:
• the different densities of forest, • small paths, large paths, roads, • marches,
• areas of forestry work,
as well as other symbols which are related to the symbols above. The international orien-teering federation has tabulated the running speed corresponding to different surfaces and veg-etation in table5.1. As we can see from the table the given percentages of how fast it is possible to run differs a lot, so to be able to calculate an optimal route choice there must be much higher accuracy in the estimate of the running speed.
No. Percentage Examples Approx. speed min/km
1 > 100% Lawns, paved areas, paths < 4
2 80 ° 100% Rough open land, forest < 5
3 60 ° 80% Stony ground, undergrowth, dense vegetation 5 ° 6 : 40
4 20 ° 60% Very stony ground, undergrowth, dense vegetation 6 : 40 ° 20
5 < 20% Extremely stony ground, very dens vegetation > 20
Table 5.1: ISOM 2017 runnability table.
anisotropy in the model. By introducing more anisotropy it is natural to implement algorithms for the more general type of anisotropic travel time equation. The fast sweeping method will for high anisotropic problems require more iterations than the stated 2niterations for problems in n dimensions in order to converge. An one-pass algorithm as the ordered upwind method would be a natural suggestion as an alternative to the fast sweeping method. Another argument for choosing an one-pass algorithm is the increasing of discretization points in the domain - when including paths and other thin symbols from an orienteering map we need high enough resolu-tion of the raster image we use to construct the speed funcresolu-tion F . An one-pass algorithm as the
ordered upwind method will restrict the computational complexity to O(N) for N2
Chapter 6
Conclusion
Level set methods appears in many applications and one of them is optimal path planning. In orienteering the athletes faces the challenges of choosing the optimal route choice on a given leg during an orienteering course, i.e. the optimal path between the control points. Today there are no tools for generating the optimal route choice on a given leg, which in the future might be an important tool for preparation before orienteering races. By introducing the problem of finding the optimal route choice in a level set framework we have different choices of algorithms for solving the resulting travel time equation.
In this thesis the isotropic travel time equation, namely the eikonal equation, was solved with the fast marching method and an anisotropic travel time function was solved with the fast sweeping method. The latter algorithm solves the problem of generating the optimal route choice on a orienteering course where the speed function F is based on a digital elevation model from LiDAR data and impassable symbols imported from an existing orienteering map in OCAD. Hence the aim of the thesis; to develop an analysis tool based on solving the travel time equation on a leg in an orienteering course, is fulfilled.
For future work on the same topic the author suggests to implement the ordered upwind method on the same problem as this algorithm have advantages in solving general anisotropic problems compared to the algorithms implemented in this thesis.
[1] F. Arnet. Arithmetical route analysis with examples of the long final courses of the world orienteering championships 2003 in switzerland and 2005 in japan. Scientific Journal of
Orienteering, Volume 17:3–20, 2009.
[2] L. C. Evans. Partial differential equations. In C. for Pure and A. Mathematics, editors, Partial
Differential Equations, volume 3A, 3B. Berkeley Mathematics Lecture Notes Series,
Univer-sity of California, Berkeley, CA, 1994.
[3] M. Hoset. Utvikle metoder for å automatisk finne beste veivalg i sprintorientering.
Høgskolen i Telemark, Prosjektoppgave i GIS, 2015.
[4] R. J. LeVeque. Numerical Methods for Conservation Laws. •, Birkhauser, Basel, 1992. [5] K.-A. Lie. an Introduction to Reservoir Simulation Using MATLAB. SINTEF ICT,
Departe-ment of Applied MAthematics, Oslo, Norway, 2016.
[6] J. R. Looker. Globally Optimal Path Planning with Anisotropic Running Costs. DSTO De-fence Science and Technology Organization, Victoria, Australia, 2013.
[7] L. C. E. M. G. Crandall and P.-L. Lions. Viscosity solutions of hamilton-jacobi equations.
Tran. AMS, pages 1–43, 1983.
[8] L. C. E. M. G. Crandall and P.-L. Lions. Some properties of viscosity solutions of hamilton-jacobi equations. Tran. AMS, pages 487–502, 1984.
[9] MATLAB. version 9.3.0 (R2017b). The MathWorks Inc., Natick, Massachusetts, 2017. [10] B. O. Myrvold. Is it possible to find a best route - a look at accuracy and significance in route
choice comparison. Scientific Journal of Orienteering, 12(1), 1996.
[11] E. Rouy and A. Tourin. A viscosity solutions approach to shape-from-shading. SIAM J.
Num. Anal, pages 867–884, 29, 3, 1992.
