Optimal Road Problem
by
Daniel Chernowitz
A thesis submitted in partial fulfillment for the degree of Bachelor of Mathematics
in the
Faculty of Mathematics and Natural Sciences Department or School Name
February 2014
Abstract
Faculty of Mathematics and Natural Sciences Department or School Name
Bachelor of Mathematics
by Daniel Chernowitz
We consider the problem of paving a path through a mountainous terrain. An optimal path is one that balances the requirement of a route to be both as short as possible and as flat as possible. A new form of Riemannian metric, ˜S, is proposed to weigh the length of a path, such as to penalize vertical transport more than horizontal, reflecting the cost and difficulty of traffic to change altitude. A geodesic γ of this metric is then a candidate for an optimal path. A number of numerical solution techniques to calculate geodesics are explored, most notably a functional iteration using Christoffel’s symbols of the geodesic equation, and geodesic evolution to minimize geodesic curvature. Finally, statistical analyses are done on the resulting geodesics of this weighted metric through a simulated hilly landscape. The metric performs well in balancing horizontal versus vertical travel economically, however it does not prevent a high maximum steepness over the course of the path.
Without the long, windy, sleep-deprived night busses from Cusco to Ica, I wouldn’t have anything to improve upon. For this reason, I thank the infrastructure of Peru. Besides this, my gratitude to professor Waalkens who has agreed to be the second corrector, for the subjects you have taught me so far. Finally, I thank my preceptor, professor Van der Schaft. A number of times, a lucid hint dissolved the walls I was trying to run through, and I have found it pleasant to work under you.
ii
Abstract i
Acknowledgements ii
1 Introduction 1
1.1 Motivation and Formulation . . . 2
1.1.1 Road Parameters and Costs . . . 2
1.1.2 Research Question . . . 4
2 Theory 5 2.1 Parameterization of the Landscape . . . 5
2.2 Weighted Path Length . . . 7
2.3 Geodesics . . . 12
3 Discrete Solution Techniques 15 3.1 Initial Estimate . . . 15
3.1.1 Dijkstra’s Algorithm . . . 18
3.2 Refinement . . . 21
3.2.1 Functional Iteration . . . 21
3.2.2 Respacing . . . 23
3.2.3 Geodesic Evolution . . . 25
3.3 Program Structure . . . 28
3.3.1 Terrain . . . 28
3.3.2 Creating Geodesics . . . 29
4 Experiment 32 4.1 Setup . . . 32
4.2 Results. . . 32
4.2.1 Distribution . . . 33
4.2.2 Discussion. . . 35
4.2.3 Comparison to Euclidean metric . . . 37
A Appendix 38 A.1 Refining Geodesics . . . 38
Bibliography 44
iii
Introduction
When travelling a winding road over mountainous terrain, the question of the efficiency arises. Should the engineers responisble have taken a different route, the road might have been shorter, or the average traffic might require less fuel to traverse it. Given certain demands on the use of the road, there must be an optimal road to take. We consider the case of a path confined to the surface of the terrain, thus disregarding bridges and tunnels, and we only require the road to pass through a source and target point on the map. Besides these conditions, the course is free.
Figure 1.1: A winding road in Peru. This segment does not need to pass any particular landmarks while navigating the terrain, so it’s course could conceivably be rerouted to
a shape that is more fuel efficient. It is subject to optimalisation. [1]
1
1.1 Motivation and Formulation
Automobiles are designed to travel horizontally. Sometimes vertical transport is in- evitable, for instance if the journey’s destination is higher than the starting point. Be- sides this, it is desirable to have as flat a course as possible, for displacing a mass against gravity constitutes performing work, which costs fuel, or resources. Even though this work is released when travelling back down, no system is 100% efficient and there will be loss. Furthermore, when vertical displacement occurs, an automobile copes best with the least possible incline. In the absurd limit, we find it cannot drive up a wall. We pro- pose that the suitability for traffic decreases as the square of the incline. Small terrain fluctuations may be appropriate, however steep hills must be discouraged when paving a way.
In this research, we will introduce a new notion of length to measure how far points are from each other on a landscape. This weighted length penalizes vertical transport.
As far as possible, we formulate it such that horizontal displacement adds to the length normally, but vertical displacement enters quadratically. Then the road between two points that minimizes this length is optimal for our purposes.
1.1.1 Road Parameters and Costs
In an attempt to quantify the demands an engineer would empose on an optimal road, we compose a formula for the total costs of placement and use. Although in reality the party responsible for construction might not need to weigh the user’s expenses as heavily as its own, here we assume the ideal situation in which they minimize total expenditure for all involved. First we set the nomenclature of quantities related to the road.
Call the length of the road L. The euclidean distance between source and target points is ∆01. When traversing in one direction, keep track of the total cumulative distance of vertical climb, Ct, and fall, Ft. Both are positive and the difference δ01:= (Ct−Ft) is the change in altitude from the source to target points. We term the unnecessary vertical travel to be the smaller of the two, U := min(Ct, Ft). Some of these measurements are illustrated in figure 1.2.
Say the road costs R in unit expense per unit length to pave. For the average vehicle, the drag against the direction of motion costs D per unit length travelled. In reality, this will vary with the instantaneous speed, and more sophisticated integrals may take braking and accelerating into account, but here we assume drivers always follow the speed limit exactly. Displacing the average vehicle a unit of length vertically costs an additional H, after correcting for the amount of potential energy released when returning downhill.
Figure 1.2: Illustration of quantities L, Ct, Ft, ∆01& δ01on a cross section of a road.
Here, U = Ft.
Finally, the road is in use for Y unit time, and the traffic consists of Q vehicles per unit time.
Definition 1.1. The total (variable) cost S of a road is a measure of how much more resources are needed to travel from source to target points than would be in the ideal situation: a straight line.
S := (L − ∆01) · (R + Y QD) + Y QU H (1.1)
The amounts Y Qδ01H and ∆01(R + Y QD) are treated as an unavoidable expenses, and the quality of the road should only be judged on the amount of excess expense. When examining a particular course, we need only know the values H, R & D to a common factor, and the total number of vehicles expected to travel the road during its lifetime, Y Q.
As rough estimates, we will use the following values.
Quantity Magnitude
R 300
D 1 · 10−4
H 2 · 10−4
Y Q 2 · 107
These are not substantiated here, as this is not a civil engineering research. The magni- tudes are derived loosely from approximate construction and gas prices, vehicle size and performances of current day, and the first three entries may be interpreted as costs in USD per meter. More appropriate values may be chosen for any particular study, these serve merely to exemplify the ruling later in this research.
As a final measure, an interesting statistic is the maximum steepness encountered. If it is too great, no matter the quality of the rest, this bottleneck will prohibit some traffic from traversing the road. This may be for vehicles lacking the traction or power to ascend or descend safely. The maximum steepness, max ∠, is defined as the maximum ratio of vertical to horizontal displacement.
1.1.2 Research Question
The research question is threefold, we state it as follows:
• Geometrically, what mathematical conditions should be satisfied by an optimal road?
• What solution techniques may be used to calculate the form of such a road?
• Given a collection of roads constructed with different parameters, how does their quality for traffic depend on these parameters and which values minimize the ex- pression in equation 1.1?
These questions will be answered in order in chapters 2, 3 and 4, respectively.
In the following, we will use the mathematical term path when referring to the equivalent of a road in an abstract landscape.
Theory
2.1 Parameterization of the Landscape
An arbitrary 3D landscape may be seen as a 2-dimensional manifold in 3-dimensional euclidean space, using cartesian coordinates. We define U :
Definition 2.1.
U (u, v) = [x, y, z] = [u, v, V (u, v)] (2.1)
U (u, v) is a 3-valued function of 2 variables. Its first and second partial derivatives are required to exist. The boundaries of the independent plane are arbitrary, but we will take V and thus U to be defined on an open connected set in R2. This construction is referred to as a Monge Patch. [2]
Definition 2.2. The first fundamental form of the landscape U , or Riemannian metric [3], is given by the following functions. [4]
E :=
∂U
∂u
2
=
1, 0,∂V
∂u
2
= ∂V
∂u
2
+ 1 (2.2a)
F := ∂U
∂u ·∂U
∂v =
1, 0,∂V
∂u
·
0, 1,∂V
∂v
= ∂V
∂u
∂V
∂v (2.2b)
G :=
∂U
∂v
2
=
0, 1,∂V
∂v
2
= ∂V
∂v
2
+ 1 (2.2c)
Arranged in a 2 × 2 matrix, because this is a 2-dimensional manifold, we obtain:
S :=
"
E F
F G
#
(2.3) 5
This defines a Euclidean metric, or distance function, on the (u, v)-plane, corresponding to the augmented space between two points due to variation in height. At each point p ∈ U , S may be used to construct an inner product between two (row) vectors l, k in the tangent space at p, TpU . [5]
hl|ki = lSkT (2.4)
In particular, this naturally defines a norm on the tangent space, and thus a vector length.
|| l || =p
lSlT (2.5)
Definition 2.3. We would now like to define the parameterization of a one-dimensional path γ in the (u, v)-plane.
γ(t) = (u(t), v(t))
t ∈ [a, b]
γ(a) = (x0, y0) γ(b) = (x1, y1)
(2.6)
We will also call (x0, y0) the source and (x1, y1) the target of the path.
This implicitely denotes a path on the 3-dimensional manifold given by the composition (U ◦ γ).
The conventional length of the path (U ◦ γ), using the metric derived in equations2.2, is the integral of the length of the velocity vector. From equation 2.5 then follows an expression for the length: [4]
L(γ) = Z b
a
E ˙u2+ 2F ˙u ˙v + G ˙v212 dt =
Z b a
˙ γS ˙γT12
dt (2.7)
V =˙ ∂V
∂uu +˙ ∂V
∂v ˙v (2.8)
An obvious result from the Pythagorean theorem is then the equality:
L(γ) = Z b
a
˙ γS ˙γT12
dt = Z b
a
˙V2+ ˙u2+ ˙v2
1
2 dt (2.9)
2.2 Weighted Path Length
Before we continue with optimization, we would like to weigh path length to penalize vertical travel more than horizontal travel. Then a path minimizing length, for our purposes called a geodesics, will avoid going up and down when possible. The formulation of geodesics will be made precise in the next section.
We would like to devise a new potential function W (u, v) and corresponding manifold UW = (u, v, W (u, v)) on which the path length of (UW ◦ γ) is the conventional path length, and equals the weighted length of γ in the landscape defined by V (u, v). If we find this W , we have reduced the question to the search for geodesics on UW.
Definition 2.4. The weighted length over the landscape of UV, subscript for clarity, is given by:
LV(γ) = Z b
a
( (∂V
∂u
2
+ 1) ˙u2+ 2(∂V
∂u
∂V
∂v) ˙u ˙v + (∂V
∂v
2
+ 1) ˙v2
12
+ α( ˙V , u, v) )
dt (2.10) It is implicit that ( ˙u, ˙v) = ( ˙u(t), ˙v(t)) = ˙γ(t). There is a linear added weight term α( ˙V , u, v), which should depend on at least the speed of vertical transport ˙V , but may also involve the position on the landscape, which is assumed known.
V is a potential function of cartesian coordinates (u, v) ∈ R2, so we may use the gradient identity
∇V = ∂V
∂u,∂V
∂v
(2.11)
Then
LV(γ) = Z b
a
n
(∇V · ˙γ)2+ || ˙γ||212
+ α( ˙V , u, v)o
dt (2.12)
Now we impose:
LV(γ) = LW(γ) = Z b
a
γS˙ Wγ˙T12 dt =
Z b a
(∇W · ˙γ)2+ || ˙γ||212
dt (2.13)
And we would like ∇W to be an analytic function of ∇V . The corresponding α will arise from the difference between UW and UV.
Theorem 2.5. UW can only take the form of a rescaling of the vertical axis in UV, i.e.
W = aV for some a ∈ R
Proof. For W to be a well defined scalar function, the vector field ∇W , its gradient, must be conservative. This is a consequence of the gradient theorem. The criterium for a vector field (X(u, v), Y (u, v)) to be conservative is equivalent to requiring its scalar curl to vanish everywhere. The scalar curl is defined as follows:
∇ × (X, Y ) = ∂X
∂v −∂Y
∂u (2.14)
Conversely, it is trivial to show that a non-conservative vector field cannot be the gradient of a potential function. On a simply-connected region, the conservative attribute is sufficient to guarantee the existence of such a W . [6] Thus we are looking for operations on ∇V that preserve the conservative feature of this vector field. Call such an operation W = (W1, W2)(x1, x2), so W = ∇W and (x1, x2) = ∇V are both conservative vector fields.
0 = ∂W2
∂u −∂W1
∂v = ∂W2
∂x1
∂x1
∂u +∂W2
∂x2
∂x2
∂u −∂W1
∂x1
∂x1
∂v −∂W1
∂x2
∂x2
∂v (2.15)
Each instance of the function W is evaluated at the point (x1, x2) = ∇V . Because of the following observations
∂x1
∂u = ∂2V
∂u2 (2.16a)
∂x2
∂v = ∂2V
∂v2 (2.16b)
∂x1
∂v = ∂x2
∂u = ∂2V
∂u∂v (2.16c)
Equation 2.15 must hold for all possible functions V (u, v). Furthermore, W is an ana- lytic vector function of two variables. It’s partial derivatives cannot be derivative opera- tors themselves. The terms ∂W∂x2
1
∂x1
∂u and ∂W∂x2
2
∂x2
∂v cannot cancel, although the remaining two can. This leads us to demand:
∂W2
∂x1
= ∂W1
∂x2
= 0 (2.17)
So
W = (W1(x1), W2(x2)) (2.18)
Substituting from equation2.16c, equation2.15 reduces to
∂
∂x1
W1(x1) = ∂
∂x2
W2(x2) = const. (2.19)
As W1 and W2 are functions of different variables, the equality can only hold if the derivatives both equal some q ∈ R. This means the only transformation W we can make is a linear one: for instance stretching the vertical axis.
∇W = W(∇V ) = q · ∇V (2.20a)
W (u, v) = q · V (u, v) (2.20b)
For some 1 6 q ∈ R. We choose q larger than unity to penalize vertical travel, not reward it.
For W defined as such, we ask what the form of α will be. Without loss of generality, we temporarily rescale the time variable to t∗ so the speed of γ(t∗) in the (u, v)-plane is unity, thus ˙u2+ ˙v2= 1. Referring back to equation 2.9, the corresponding cost function α takes the following form:
α( ˙V , q) = γS˙ Wγ˙T12
− ˙γSVγ˙T12
= q
(q ˙V )2+ 1 − q
( ˙V )2+ 1 (2.21)
In the limits ˙V = +∞ and ˙V = −∞, α( ˙V ) will go to q
(q ˙V )2+ 1 − q
( ˙V )2+ 1 ≈ q
(q ˙V )2− q
( ˙V )2 = |a ˙V | − | ˙V | = (q − 1)| ˙V | because q > 1. The form of the penalty is illustrated in figure 2.1.
Though searching for distance-minimizing paths on this rescaled landscape would pro- duce paths biased against vertical transport, it is not very innovative from a mathemat- ical point of view. This amounts to finding geodesics on a new euclidean surface, which has been done in many applications before.
In light of the result of theorem 2.5, we abandon the search of a landscape function on which normal path length of γ equals the weighted path length with quadratically penalized vertical transport on U . However, an alternative Riemannian metric may be proposed to calculate this weighted path length. For this metric, there may be no well defined landscape that produces it, as in equation 2.2, however, nothing fundamental prohibits us from calculating Christoffel symbols and finding geodesic curves.
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0
1 2 3 4
V˙
α
α(q = 3)
|2 ˙V |
Figure 2.1: Cost function α as a function of ˙V . We have taken q = 3, and far from the origin α approaches 2| ˙V | asymptotically.
From the formulation of S in definition2.2, we observe that it may be represented as:
S = T (V, u, v) + I2×2=
"
Vu2 VuVv
VuVv Vv2
# +
"
1 0 0 1
#
(2.22)
Using the shorthand Vu := ∂V∂u and Vv := ∂V∂v
Here T is dependent of the shape of the landscape V , observe that ˙γT ˙γT = ˙V2. The identity component I2×2 stems from the fact that u & v are cartesian orthogonal co- ordinates that always contribute to the travelling speed by the Pythagorean theorem:
˙
γI ˙γT = ˙u2+ ˙v2.
Definition 2.6. To weigh paths against steep vertical travel, we propose the alternative shape operator or metric: ˜Sq = qT2+ I.
Expanding:
S˜q= q
"
Vu4+ Vu2Vv2 Vu3Vv+ Vv3Vu Vu3Vv+ Vv3Vu Vv4+ Vu2Vv2
# +
"
1 0 0 1
#
= q||∇V ||2T + I (2.23)
The motivation behind this definition is that the steepness of the landscape enters quadratically. Additionally, we may weigh this dependency through the variable a.
The scalar function ||∇V || is the absolute steepest rate of ascent/descent at each point of U . Note we will often drop the subscript q to streamline equations.
We proceed by exploring expressions of the path length of γ, using this metric.
LV(γ) = Z b
a
˙
γ ˜S ˙γT12 dt =
Z b a
q||∇V ||2γT ˙˙ γT + ˙γ · ˙γT12 dt
= Z b
a
q||∇V ||2V˙2+ ˙u2+ ˙v212 dt
(2.24)
The second equality is due to the distributive property of matrix multiplication.
We now wish to describe the behaviour of the cost function, α( ˙V , u, v). As before, we assume the speed of γ to be unit.
α( ˙V , q, u, v) =
γ ˜˙S ˙γT
12
− ˙γS ˙γT12
=
q||∇V ||2V˙2+ 1
12
− ˙V2+ 1
12
(2.25)
When travelling along the direction of steepest ascent/descent with unit speed, ||∇V || =
| ˙V |. Equation2.25reduces to α = q
q ˙V4+ 1 −p ˙V2+ 1. The shape of the cost function is illustrated in figure2.2.
−6 −4 −2 0 2 4 6
0 10 20 30
V˙
α
q = 1 q = 2
Figure 2.2: Cost function α as a function of ˙V , for a path along the steepest ascen- t/descent. We have taken parameters q = 1 and q = 2.
We see the function is flattened around the origin but adopts a quadratic trend for higher vertical speeds. This quadratic trend diminishes if the travel is not along the steepest ascent. A quirk of this metric is that the minimum cost occurs at a nonzero value of ˙V , although the difference |αmin− α( ˙V =0)| is small by physical standards.
If the angle between the direction of steepest ascent, and the direction of unit speed travel is θ, then the magnitude of the instantaneous vertical speed is given by | ˙V | =
cos θ·||∇V ||. Inserting this into equation2.25allows us to visualize the case of movement in all directions. See figure2.3.
−π2
0
π 2
−4 −2 0 2 4
0 20 40
θ
||∇V ||
α
q = 2
Figure 2.3: Cost α as a function of the local maximum steepness ||∇V || and the angle θ between the direction of γ and the direction of maximum steepness, for q = 2.
2.3 Geodesics
As explained in chapter 1, we will be searching for paths to connect two points that in some ways minimize the path length between them. We see the path length in equation2.24is an integral over an expression involving an undetermined function. This expression is a functional, and an extreme value of the integral can be obtained through the methods of calculus of variations. In particular, the Euler-Lagrange equations solve this problem. [? ] However, in this research, we will use the theory of geodesics, which are locally distance minimizing. Formally, that means the shortest path over a manifold connecting two points lies along a geodesic curve. It can be shown that these are the same solutions as obtained by the Euler-Lagrange equations. [3]
Geodesic curves have the property that their speed, though not necessarily velocity vector, is constant and unit. Any curve can however be reparameterized to meet this requirement.
d
dt||(U ◦ γ(t))|| = 1 (2.26)
More importantly, a geodesic’s acceleration is always normal to the surface of the man- ifold. For some C ∈ R:
d2
dt2(U ◦ γ(t)) = CN = C∗ ∂V
∂u,∂V
∂v, −1
(2.27)
Normalization of surface normal N is absorbed into the second constant C∗. This condition together with that of equation2.26are tantamount to the following differential equations:
¨
u + Γuuuu˙2+ 2Γuvu u ˙v + Γ˙ uvv˙v2 = 0 (2.28a)
¨
v + Γvuuu˙2+ 2Γuvv u ˙v + Γ˙ vvv˙v2 = 0 (2.28b)
The functions Γ••• are called the Christoffel symbols of the manifold. They can be calculated by the matrix equations: [7]
"
Γuuu Γvuu
#
=
"
E F
F G
#−1" 1
2Eu
Fu−12Ev
#
(2.29a)
"
Γuuv Γvuv
#
=
"
E F
F G
#−1"1
2Ev
1 2Gu
#
(2.29b)
"
Γuvv Γvvv
#
=
"
E F
F G
#−1"
Fv−12Gu
1 2Gv
#
(2.29c)
In the above, the functions E, F and G are as defined in equations 2.2. Also, we have used the shorthand ∂E∂u = Eu and similar for the other components. Finally, [ • ]−1 denotes the matrix inverse, given by:
"
E F
F G
#−1
= 1
EG − F2
"
G −F
−F E
#
(2.30) Definition 2.7. A γ(t) that solves equations2.28is a geodesic, which locally minimizes distance with respect to the Riemannian metric S.
Putting things together, we will now construct the geodesic equation for our optimal path, as measured by the weighted metric of definition2.6. Equations 2.29produce the Christoffel symbols. To this end we insert the identities:
S˜−1 = 1 1 + q||∇V ||4
"
q||∇V ||2Vv2+ 1 −q||∇V ||2VuVv
−q||∇V ||2VuVv q||∇V ||2Vu2+ 1
#
(2.31)
Which follows from equation2.30, and:
Eu = 2q (2Vu3+ Vv2Vu)Vuu+ Vu2VvVuv
(2.32a) Ev = 2q (2Vu3+ Vv2Vu)Vuv+ Vu2VvVvv
(2.32b)
Fu = q (3Vu2Vv+ Vv3)Vuu+ (3Vv2Vu+ Vu3)Vuv
(2.32c)
Fv = q (3Vu2Vv+ Vv3)Vuv+ (3Vv2Vu+ Vu3)Vvv
(2.32d)
Gu = 2q (2Vv3+ Vu2Vv)Vuv+ Vv2VuVuu
(2.32e)
Gv = 2q (2Vv3+ Vu2Vv)Vvv+ Vv2VuVuv
(2.32f)
These are found by direct calculation. Finally, we define the matrices Γu and Γv:
Γu=
"
Γuuu Γuuv Γuuv Γuvv
#
, Γv =
"
Γvuu Γvuv Γvuv Γvvv
#
(2.33)
Using this substition, the geodesic equations of2.28 take the form
¨
γ = − ˙γΓuγ˙T, ˙γΓvγ˙T
(2.34)
Another classification of a geodesic is that it is a curve for which the geodesic curvature κg vanishes everywhere. Classically, we may view κg as a measure of the acceleration of γ(t) in the tangent plane to the surface. With our altered metric, the notion of the tangent plane is not defined here, but this curvature may nonetheless be formulated in terms of the metric and Christoffel symbols. [2] We shall find κg useful later on.
κg :=
v u u u t
det( ˜S) h
˙ γ ˜S ˙γT
i3 · ˙γ( ˙vΓu− ˙uΓv) ˙γT + ¨u ˙v − ¨v ˙u
(2.35)
Discrete Solution Techniques
We now move on in an attempt to solve equation 2.34. As there will not be an analytic solution for the general case of an arbitrary landscape defined by V (u, v), we turn our attention to approximation techniques. If we discretize the environment, restricting γ to a succession of straight line segments passing through a finite set of points in the (u, v) plane, then it is possible to determine the optimal path under these conditions. In this chapter, we make this notion concrete.
3.1 Initial Estimate
In graph theory, consider a graph with weighted edges. There are a finite number of ways to travel from a starting vertex to a goal vertex, granted we forbid self-intersection.
This may also be done in our case by creating a discrete graph to describe the contin- uous landscape, to a good degree of approximation. As a first step in discretizing the environment, we define a mesh, analogous to that in [5].
Definition 3.1. Let U (u, v) be defined at least inside the rectangular domain (u, v) ∈ [uB, uC] × [vB, vC]. Then for n ∈ N>1, the mesh Mn is a set of (n + 1)2 evenly spaced points, of this domain, positioned on a rectangular grid. We assume the boundaries of the domain of U are sufficiently far from the target and source that expanding the domain further does not facilitate a better path.
Mn= (ui, vj) =
uB+(uC− uB)
n i, vB+(vC− vB)
n j
i, j ∈ {0, 1, ..., n} (3.1)
15
This mesh represents a graph, if we treat each grid point (ui, vj) as a vertex, connected by edges to each of its direct or diagonal neighbors. The terms grid point and vertex will be synonymous in the following.
Definition 3.2. The induced graph on the domain (u, v) ∈ [uB, uC] × [vB, vC] is the pair of sets (Mn, Bn) in which Bn consists of subsets with two elements of Mn that are connected by an edge.
Bn= ((ui, vj), (ui+k, vj+l))
(ui, vj), (ui+k, vj+l) ∈ Mn
k, l{−1, 0, 1} (3.2)
An example of such a graph is shown in figure3.1.
Figure 3.1: An example of a mesh with n=6. The boundaries of the domain are noted, as well as a point(ui, vj). Also, in the bottom right, the edges that originate
from any element of M6 not on the boundary are drawn.
Definition 3.3. A discrete path γ∗(t) = (u∗(t), v∗(t)), t ∈ [a, b] is the continuous, piece- wise differentiable function obtained by subsequently connecting neighboring elements of Mn by straight lines, without self-intersections. Furthermore, γ∗(a) = (x0, y0) &
γ∗(b) = (x1, y1). This places a slight constraint on the possible geodesics the path can mimic, but with a sensible choice of uC, uB, n etc. this should not be an obstacle. For an illustration, see figure3.2.
A note on notation, because the time-parameterization proves unwieldy. If γ∗ passes through r grid-points, then γ∗h ∈ Mn is hth vertex of the path, h ∈ {1, 2, ..., r}. In particular, γ∗(a) = γ∗1 and γ∗(b) = γ∗r
Figure 3.2: A discrete path γ∗(t) on M6. Certain mesh-points on the path are signified.
Finally, we wish for this graph to approximate traversing the landscape. We construct a weight function on Bn, such that the length of the edges mimics the metric ˜S, as formulated in equation 2.24.
Definition 3.4. The weight function ωq : Bn → R is constructed as a best discrete approximation to the continuous ˜S. It is given by the following formula:
ωa: ((ui, vj), (ui+k, vj+l)) 7→
q
(∆u, ∆v) ˜S(q, ¯u, ¯v) (∆u, ∆v)T (3.3)
Here, ∆u and ∆v are the horizontal step length over this edge in the x and y direction, respectively. The arguments of ˜S signify that the metric is evaluated at the midpoint between (ui, vj) and (ui+k, vj+l) in the (u, v)-plane.
For a path traced over an ordered set of points on the graph, this weight function presents a sensible discrete appoximation of the weighted path length. L∗V is defined to equal the sum of the weights of the traversed edges.
3.1.1 Dijkstra’s Algorithm
Dijkstra’s algorithm is an efficient tool to find the shortest path through a grid. It constructs a tree, so that for each vertex in Mn the weighted length of the shortest connection to source vertex γ∗1 = (x0, y0) is tabulated. When this cartography is com- plete, say the shortest path to the target γ∗r = (x1, y1) has length L. Backtracking from the target vertex, the direct neighbor with the shortest path to the source is ob- viously on the path of length L. Iteratively, continuing to hop to the neighbor with the shortest tabulated path length will eventually connect the target to the source, and this sequence defines the optimal path γ∗. [8] The process of tabulating the smallest distance to each vertex, modified for a graph induced on a mesh, is described in the box denoted Algorithm 1.
Result: Tabulating the shortest distance from the source vertex to each other vertex For mesh Mn, construct two (n + 1) × (n + 1) matrices, D with all entries equal to ∞ and E , with all entries zero. D contains the tentative shortest distance to each vertex, and E keeps track of which vertices have been visited so far. Indices begin numbering at zero, so D0,0 is the top left element of D, and in general elements are called by Di,j for i, j ∈ {0, 1, ..., n}, analogously for E . It is understood that Di,j and Ei,j correspond to the vertex ui, uj of Mn;
Set the current element (I, J ) to correspond to the source vertex, i.e. (uI, vJ) = (x0, y0). Each element will only be current once. Afterwards, the element or vertex is said to be visited ;
Additionally, set DI,J := 0, the distance from the source to itself is 0;
while The target vertex is unvisited: Ei,j 6= ∞, where (ui, vj) = (x1, y1) do
For all neighbors of the current element, check whether the path coming directly from the vertex of the current element is shorter than the distance already tabulated to that neighbor: for k, l ∈ {−1, 0, 1} do
if The neighbor is on the grid, (I + k), (J + l) ∈ {0, 1, ..., n} then δd:= min{DI+k,J +l, DI,J+ ωa((uI, vJ), (uI+k, vJ +l))};
If applicable, replace the tentative distance at the neighbor by the new, shorter distance. DI+k,J +l:= δd;
end end
Mark the current element as visited. EI,J := ∞;
Reset the current element (I, J ) as that corresponding to the unvisited vertex with the smallest tentative distance, such that (D + E )I,J := min
i,j {Ei,j+ Di,j}.;
end
Algorithm 1: Dijkstra’s algorithm on a discretized landscape
One may note that not all vertices have to be consideren by this algorithm. If all unvisited vertices have a larger tentative distance to the source than the target, a path through them cannot be shorter than one through the collection of visited vertices.
We will not go into detail proving the working of dijkstra’s algorithm, as it is a well known and frequently used tool.[8] We observe that in this implementation, there will always be multiple connecting paths amongst which to minimize weighted length. Take, for instance, one commencing with solely horizontal travel to align the path with the u-coordinate of the target, and then vertical travel to align the v-coordinate.
We do not know a priori the resulting discrete path γ∗of Dijkstra’s algorithm. Thus it is impossible to judge in the general case how well U ◦ γ∗ will approximate a geodesic U ◦ γ on the same surface. We can however construct a different discrete path to stay close to a smooth curve. The question arises whether a discrete path’s weighted length can be comparable to that of an arbitrary, given geodesic. We claim that for a sufficiently fine mesh, this is the case, and we will put bounds on the ratio between the path lengths.
Then, because both the geometrical as well as the combinatorial approach produce an optimal path, we know that their respective results will not lie far apart.
Theorem 3.5. In the limit n → ∞, a discrete path γ∗over the mesh Mnsharing source and target points with a smooth curve γ, can be made to satisfy LV(γ∗) 6 CLV(γ), for a suitable constant C ∈ R.
Proof. Consider a rectangular patch P of (m + 1) × (m + 1) grid points of the mesh that γ = (u(t), v(t)) passes through, m ∈ N>1. This patch has dimensions mn([uC − uB] × [vC − vB]). For suffiently large n, and thus sufficiently fine mesh, we may make the approximation that inside P , γ is a straight line, and ˜S is constant. Furthermore, without loss of generality, we confine our attention to the case where 0 6 dvdu = ˙v/ ˙u 6 1.
Because the edges on the grid are symmetrical in all diagonal and right-angle directions, any scenario in which the slope of γ lies outside the interval [0, 1], in another octant, is analogous to this one. We then construct a test path γ∗ as follows: On the left edge of the patch P , take the first point of γ∗ on the vertex closest to where γ enters P . Then connect this point to a vertex to the right along the u-direction, ∆v = 0 or diagonally up to the right, ∆v 6= 0, such that within that column of gridpoints γ∗ and γ are closest together. Iterate this procedure m times, until the patch is traversed. For an illustration, see figure3.3.
Within P , because the metric ˜S is taken constant, we may rearrange the diagonal and straight segments of γ∗to be grouped together. The number of straight edges is denoted l, the number of diagonal edges is k and as an approximation, we will assume that inside P , γ also travels the length of k edges in the v-direction. The weighted length LV is
Figure 3.3: The smooth curve γ is locally assumed linear and runs through a small patch P of the mesh. The construction of a test path γ∗ such that it remains close to
γ is illustrated.
defined in equation2.24. The ratio of the discrete to smooth lengths in the patch is then given by:
LV(γ∗) LV(γ)
P
= Rl
0
q
(1, 0) ˜S(1, 0)Tdt +Rk
0
q
(1, 1) ˜S(1, 1)Tdt R(l+k)
0
q
(1,k+lk ) ˜S(1,k+lk )Tdt
(3.4)
Here, we have parameterized both paths by their u-coordinate and have used := mn(uC− uB), the ∆u between grid points. The elements of ˜S may be called by ˜S =
"
s11 s12 s12 s22
# . Substituting the auxiliary variable λ := (l/k) and Σ ˜S := (s11+ 2s12+ s22), we obtain:
LV(γ∗) LV(γ)
P
= λ√
s11+p Σ ˜S
ps11(λ + 1)2+ 2(λ + 1)s12+ s22 (3.5) In the limiting cases λ = 0 and λ = ∞, this expression equals unity, as expected. In between it exhibits an extreme value, which defines C1, itself a function of the metric.
C1( ˜S) := max
λ∈R>0
LV(γ∗) LV(γ)
P
= 2s11(Σ ˜S) − 2(s11+ s12)p s11Σ ˜S det( ˜S)
!12
(3.6)
We observe that det( ˜S) = 1 + a2||∇V ||4Vv2Vu2 > 1, so this expression is well defined and finite. We know from the triangle inequality that this ratio can never drop under unity, so equation 3.6 must be a maximum. The metric ˜S need not exhibit a symmetry in 8 directions. Analogously to C1, we define Ci, i ∈ {2, 3, ..., 8} to be the ratio of discrete to smooth weighted length, for the cases where the slope of γ|i lies in the other 7 octants.
Finally, to complete the proof, we define C to be the maximum value attained by Ci( ˜S) over all i and the domain of U . Then for the whole smooth path γ and it’s approximating path discrete counterpart γ∗, no matter the direction:
C = max
i max
P ⊂Mn
LV(γ∗) LV(γ|i) P
⇔ LV(γ∗) 6 CLV(γ) (3.7)
3.2 Refinement
The result of Dijkstra’s Algorithm is sequence of straight paths, none of which are necessarily geodesic segments. This forms a first estimate to a geodesic, though, and the methods described below may be used to refine it. From here, we abandon the notion that the points of γ∗lie on the grid, as these they may be any ordered collection points in the (u, v)-plane. The refinement methods will perturb the position of these points, and their resulting positions may take any values in R2, inside or outside [uC, uB] × [vC, vB].
We assume U (u, v) is still well defined on γ∗. Also, the discrete path length L∗V is generalized to weigh each line segment analogous to equation 3.3, with ˜S evaluated halfway between the endpoints of the segment.
3.2.1 Functional Iteration
We would like the path produced to satisfy the geodesic equation. Because γ∗ is a discrete collection of points in a plane, we must construct a discrete condition that is analogous to equation 2.34. Mimicking the technique in [5], we make the following observations. For γ∗h= (uh, vh):
˙
uh= uh+1− uh−1
∆th+ ∆th−1
+ O(∆t2) (3.8a)
˙vh= vh+1− vh−1
∆th+ ∆th−1 + O(∆t2) (3.8b)
¨
uh= 1
∆th+ ∆th−1
uh+1− uh
∆th −uh− uh−1
∆th−1
+ O(∆t2) (3.8c)
¨
vh= 1
∆th+ ∆th−1
vh+1− vh
∆th −vh− vh−1
∆th−1
+ O(∆t2) (3.8d)
In the equations above, ∆th is the time-step between γ∗h and γ∗h+1. Because a geodesic is arclength parameterized, this is also the weighted length between these points. We assume that for each h, these are of comparable magnitude, represented by ∆t.
If one can ensure that the time-steps between all consecutive points of γ∗ are sufficiently uniform and equal to ∆t, these may be divided out. We make the approximation that higher-order terms may be truncated.
˙
uh≈ uh+1− uh−1
2∆t := ∆uh
2∆t (3.9a)
˙vh≈ vh+1− vh−1
2∆t := ∆vh
2∆t (3.9b)
¨
uh≈ uh+1+ uh−1− 2uh
∆t2 (3.9c)
¨
vh≈ vh+1+ vh−1− 2vh
∆t2 (3.9d)
We define ∆uh := uh+1− uh−1 and ∆vh equivalently.
Inserting into equation 2.34, we obtain a pair of relations to satisfy for every point of γ∗, one for the u and v coordinate, respectively.
uh+1+ uh−1− 2uh
∆t2 = − ∆uh 2∆t,∆vh
2∆t
Γu ∆uh 2∆t,∆vh
2∆t
T
(3.10a) vh+1+ vh−1− 2vh
∆t2 = − ∆uh 2∆t,∆vh
2∆t
Γv ∆uh 2∆t,∆vh
2∆t
T
(3.10b)
Implicitly, the matrices γ• are evaluated at γ∗h. We eliminate the time steps and formu- late expressions for uh and vh, in terms of their neighbors. For additional brevity, set the vector ∆γ∗h:= (∆uh, ∆vh).
uh = uh+1+ uh−1
2 +1
8(∆γ∗h)Γu(∆γ∗h)T (3.11a) vh = vh+1+ vh−1
2 +1
8(∆γ∗h)Γv(∆γ∗h)T (3.11b)
ξ∗h := (uh, vh) (3.11c)
N.B. this result has a different form in [5], which the author believes to be erroneous.
We employ the equations3.11 in a method of refinement.
Definition 3.6. Define the method of functional iteration on the geodesic. For γ∗ consisting of r ordered points, construct the function F : R2r−4→ R2r−4, for these sets are the domain of the inner, variable points of γ∗.
F : (uh−1, uh+1, vh−1, vh+1) 7→ ξ∗h | h ∈ {2, . . . , r − 1} (3.12)
The points (u1, v1) = (x0, y0), (ur, vr) = (x1, y1) are stationary and are not explicitely counted as arguments of F . Iteratively applying F to the discrete path will be referred to as functional iteration.
A discretized, true geodesic ξ∗ with uniform time-steps forms a fixed point of this func- tion: F (ξ∗) = ξ∗, as the geodesic equations 3.11are satisfied. Empirically, this method has the effect of smoothing γ∗and, importantly, reducing the magnitude of L∗V(γ∗) with each iteration, suggesting the path converges on a true geodesic.
3.2.2 Respacing
The stability of functional iteration above proves to rely crucially on a uniform time- step between subsequent points of the discrete path. For this reason, we develope an algorithm for spacing to interpolate a discrete path with evenly spaced points.
Definition 3.7. A discrete path γ∗ consisting of r points has a spacing of δs, if the following is satisfied.
q
(uh+1− uh, vh+1− vh) ˜S(uh+1− uh, vh+1− vh)T = δs |∀h ∈ {1, . . . , r − 1} (3.13)
To facilitate the construction of this spacing, ˜S is evaluated at (uh, vh). See algorithm 2.
Result: Interpolating a discrete path γ∗ with spacing δs, resulting in discrete path ξ∗ δs must be significantly smaller than the minimum arc length between any two
consecutive points of γ∗;
The first points coincide, ξ∗1 := γ∗1; H := 2, the current element of ξ∗; for each point of γ∗ : h ∈ {1, . . . , r − 1} do
Choose ξ∗H on the line between γ∗h & γ∗h+1, such that it is arc length δs from ξ∗H−1. δs=
q
(ξ∗H − ξ∗H−1) · ˜S(ξH−1∗ ) · (ξ∗H − ξ∗H−1)T. This amounts to solving for the positive root of a quadratic equation with positive discriminant, we elaborate on this calculation in equation3.16;
Set the distance to the next point of γ∗ in the (u, v)-plane, R := ||γ∗h+1− γ∗h||;
Calculate ΣR:= ||ξ∗H − γ∗h||, the distance covered by ξ∗ from the current point of γ∗ toward the next;
while ΣR< R do
H := H + 1, the next point of ξ∗ becomes current;
Set θ := arctan(uvh+1−vh
h+1−uh), correcting for quadrant such that e1:= (cos θ, sin θ) is the unit vector pointing from γ∗h to γ∗h+1. Place ξ∗H in the (u,v)-plane on the line on the line spanned by γ∗h & γ∗h+1, an arc length δs from ξ∗H−1:
ξ∗H := ξ∗H−1+ δs/ q
e1· ˜S(ξ∗H−1) · eT1; ΣR:= ||ξ∗H − γ∗h||;
end
The last point placed, ξ∗H, surpassed the original path, but it’s position will be corrected in the next iteration of the for-loop.
end
Algorithm 2: Respacing algorithm for a discrete path
It will occur that subsequent points of ξ∗lie in different segments of γ∗. In general, these segments will not be parallel, so ξ∗ must turn the corner. The situation is illustrated in figure3.4.
Define ∆ξ := (γ∗h− ξ∗H−1) to be the vector component ξ∗ must still cover in the (u, v)- plane before γ∗h. Now we call ρ the euclidean distance to be covered between γ∗h and ξ∗H), in the direction of (γ∗h+1. Also, e1:= (cos θ, sin θ) is constructed as in algorithm 2.
The following must hold:
δs= q
(∆ξ + ρe1) · ˜S(ξ∗H−1) · (∆ξ + ρe1)T (3.14)
Now, we cannot factor out ρ as is possible if all points are collinear. We employ the quadratic formula with a discriminant, Dρ.
Figure 3.4: ξ∗H−1& ξ∗H are on opposite sides of the current point γ∗h
δ2s = (∆ξ + ρe1) · ˜S · (∆ξ + ρe1)T (3.15a) 0 = e1Se˜ T1ρ2+ (∆ξ ˜SeT1 + e1S(∆ξ)˜ T)ρ + ∆ξ ˜S(∆ξ)T − δ2s (3.15b) Dρ= (∆ξ ˜SeT1 + e1S(∆ξ)˜ T)2− 4e1Se˜ T1(∆ξ ˜S(∆ξ)T − δs2) (3.15c)
ρ =
−
∆ξ ˜SeT1 + e1S(∆ξ)˜ T
+pDρ
2e1Se˜ T1 (3.15d)
By construction, equation3.15 has real solutions, so Dρis positive. Since also distance is constructed to be nonnegative, we select the positive root in the final line. Note that any metric is positive definite [3], thus the denominator is positive. With the magnitude of ρ known, we simply place ξ∗H displaced by ρ along the line towards γ∗h+1.
ξ∗H := γ∗h+ ρe1 (3.16)
As a final remark, uniform spacing comes at the cost of perturbing the target point γ∗r, by at most δs. Although this may be acceptable, one must be wary not to allow the target point to drift further by applying the algorithm multiple times unchecked.
3.2.3 Geodesic Evolution
After functional iteration, we explore a second method due to [9] for improving the dis- crete approximation of a geodesic. It involves calculating the local geodesic curvature
