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 . . . 7

2.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

1_{The map used and route choices drawn in}_{Figure 1.1}_{are collected from}_{http://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 = FdT _{d 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 (x*0*, t*0), then

*vt(x*0*, t*0*) + H(@v(x*0*, t*0*), x*0) ∑ 0 (2.5)

*2. if u ° v has local minimum at a point (x*0*, t*0), then

*vt(x*0*, t*0*) + H(@v(x*0*, t*0*), x*0) ∏ 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.7*is 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 un _{i}*.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}_{+}*k*2
2 *uxx+ o(k*
3_{),} _{(2.8)}

and this equation can be formulated (by following the notation in [12]) as the following
scheme
*u _{i}n+1_{= u}n_{i}*

_{° kD}0xu_{i}n_{+}

*k*2 2

*D*

*+x°x*

_{u}n*i*, (2.9)

*where D0x*

_{u}n*i* is shorthand for*u(x+¢x,t)°u(x°¢x,t)2¢x* *and D+x°xuin*is 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 that*d x _{d t}*

*= u we identify*Equation 2.13

thus the right hand side is zero and *du _{d 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

*Ux*2*+Uy*2*+Uz*2° 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,*

*u _{i +1}n*

*° un*

_{i}¢*t* =

*g (u _{i}n,u_{i +1}n*

*) ° g(u*

_{i °1}n*,un*)

_{i}¢*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°xun _{i}* 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,

*u _{i}n+1_{= u}n_{i}*

_{° ¢t[max(0, ai})D°xun_{i}*]. (2.21)*

_{+ min(0, ai})D+xu_{i}n*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°x*

_{u}n*i* *is used. And, consequently, when a < 0 the expression*
*max(0, ai) = 0 and min(0,ai) = ai* *thus the difference operator D+x _{u}n*

*i* is used. This idea can be
used to rewrite theEquation 2.2as,

2
6
6
4
*max(D°x _{u}n*

*i j k*,0)2

*+ min(D+xuni j k*,0)2

*max(D°y*

_{u}n*i j k*,0)2

*+ min(D+yuni j k*,0)2

*max(D°z*

_{u}n*i j k*,0)2

*+ min(D+zuni j k*,0)2 3 7 7 5 1 2 =

*1*

_{F}*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 ,j*to 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 ,j*to 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 in*Figure 2.3*at a grid point ¯x*0. 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 ¯x*0*with constant rs,*

a possible travel-time function has the form,

*T ( ¯x) =*q*( ¯x ° ¯x*0)*TK*°2*( ¯x ° ¯x*0*) ° g( ¯x ° ¯x*0)*Tq,*¯ (2.24)
where,
*• K = Q*
√
*a 0*
*0 b*
!
*QT, Q = [ ¯q, ¯q*?],
*• a =*≥1_{2}≥* _{v}*1

*+*

_{u}*v*1

*d*¥¥

_{°1}

*, b = vt*and

*• g =*1

_{2}≥

*1*

_{v}*°*

_{d}*v*1

*u*¥ (gravitation term)

We now have the framework for the alternative travel time equation. Taking the gradient of

*T gives,*
*rT =* p 1
*( ¯x ° ¯x*0)*TK*°2*( ¯x ° ¯x*0)
*K*°2*( ¯x ° ¯x*0*) ° g ¯q*
= 1
*|K*°1_{( ¯x ° ¯x}_{0}_{)||}*° 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 = v**u*= 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°x _{T}n*

*i j k,°D+xTi j kn*,0)2

*+max(D°yT*

_{i j k}n*,°D+y*

_{T}n*i j k*,0)2

*+max(D°zT*

_{i j k}n*,°D+zT*,0)2 3 7 7 5 1 2 =

_{i j k}n*1*

_{F}*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. T _{i +1,j,k}, T_{i °1,j,k}, T_{i ,j +1,k}, T_{i ,j °1,k}, T_{i ,j,k+1}, T_{i ,j,k°1}*. We assume all the neighbor points are
known, so we solve equation (3.1

*) for Ti ,j,k. Assuming we have N*2grid 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

*N*2discretization 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 u _{i ,j}h*

*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,*

*[(uh _{i ,j}_{° u}h_{xmi n}*)+]2

*)+]2*

_{+ [(u}h_{i ,j}_{° u}h_{ymi n}*2*

_{= f}_{i ,j}*hh*(3.3)

*where in the x-direction uh*

_{xmi n}= min(u_{i °1,j}h*,uh*)

_{i +1,j}), in the y-direction uh_{ymi n}= min(uh_{i ,j °1},u_{i ,j +1}h*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 uh _{i ±1,j}*

*and uh*

_{i ,j ±1}followed by updating u_{i ,j}h*according to uupd ated*

_{i ,j}*, ¯*

_{= min(u}_{i ,j}old*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(uh _{xmi n},uh_{ymi n}) + fi ,jh*

_{|u}h_{xmi n}_{° u}h_{ymi n}_{| ∏ f}i ,jh*uh*

*xmi n+uhymi n*+

q
*2f*2

*i ,jh*2*°(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 2*n* *sweeps for problems in n-dimensions. Now, as we are speaking of more general cost*

functions, we might need more than 2*n*sweepings 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 x*0*, and an end point x*1, 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 x*1. It takes the value from

*the pixel which x*0*represents, and then all the other values in W (x) are given by taking the *

*dif-ference between the pixel value in x*0and 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 point1_{and }

*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 x*1.

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) = x*1*and ø¬> 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 ∞k*is the approximation ofEquation 4.1.1*at ¬ = øk.*

FromFigure 4.3we see the geodesic curve that connects the departure- and arrival
*coor-dinate, respectively x*0*and x*1. The geodesic curve does not move perfectly through the large,

*Figure 4.3: Figure showing the geodesic path between x*0*and x*1.

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*.

*G*0 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-2_{LiDAR data in Norway is distributed through}_{https://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
*x*2* _{+y}*2

_{°r}_{0}¥

_{(}

*2*

_{ax}*2*

_{+2bx y+c y}*) °1ØØØ*

_{+d}_{Ø}

_{(4.5)}

*where r*0*= 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.1x*2* _{+ x y + y}*2

_{+ 0.4)}p

*x*2

*2*

_{+ y}*+ (x + 2y)*µq

*x*2

*2*

_{+ y}_{° 0.2}∂! · ...

*...e*≥

_{p}

*x*2

*2*

_{+y}_{°0.2}¥

_{(}

*2*

_{0.1x}*2*

_{+x y+y}_{+0.4}) ·... ... µ

*e*≥

_{p}

*x*2

*2*

_{+y}_{°0.2}¥

_{(}

*2*

_{0.1x}*2*

_{+x y+y}_{+0.4}) °1∂±µ

*e*≥

_{p}

*x*2

*2*

_{+y}_{°0.2}¥

_{(}

*2*

_{0.1x}*2*

_{+x y+y}_{+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 e*_{1}_{=} *k ˆu°uexact*k1

*kuexact*k1 . 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(*_{40}*x* ) + 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 pu*from,

*• 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 2*niterations 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 N*2

**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.*