Skip to topic
|
Skip to bottom
Jump:
TModeling
TModeling Web
TModeling Web Home
Changes
Notify
Index
Search
Webs
BioGeometry
Main
TModeling
TWiki
Edit
Attach
Printable
TModeling.ShortestPathComparison
r1.1 - 23 May 2005 - 01:49 - Main.guest
topic end
Start of topic |
Skip to actions
---+ Comparing shortest path methods ---++ Overview We will discuss the problem of finding the shortest path between two points on a terrain. Terrains can be represented in many formats. Two common terrain representations are raster representations and graph representations. Raster representations store the elevation of the terrain at regularly spaced points. Graph representations store the terrain as irregularly spaced sample points in 3D space. These sample points form the vertices of the graph; the edges are created through Delaunay triangulation. One common graph representation for terrains is a TIN (triangulated interconnect network). The shortest path between two points is easy to find on a perfectly flat terrain with no obstructions; in this case, the shortest path is simply a straight line between the two points. In mountainous terrain, however, the shortest path between two points is not as easy to find. When computing paths over terrain where the elevation varies, we can consider the shortest path between two points to be the path on the surface of the terrain that has the shortest distance. To do this we sample the elevation at various points, giving us points in 3D space. We then use the distance in 3D space to these sample points as our distance. This is the simplest option. We may want to weight paths to certain points to indicate the relative difficulty of traveling between these points. We might consider going uphill harder than going downhill and assign a greater weight to an increase in eleveation than to an equal decrease in elevation. We might want to avoid steep sections, so we can use slope to place high weights on steep sections so that a shortest path will be less likely to use them. We may have routes that should not be considered, such as those that go over cliffs, and can add infinite weights to these routes. We consider the problem of finding the shortest path in 3D space that is on the surface of the terrain. The following is a brief summary of two approaches to do so. ---+++ Sack and Lanthier Sack and Lanthier [2] describe an approach to compute a shortest path on triangulated terrain. This approach relies on the fact that there are well-known algorithms, such as Dijkstra's algorithm, that can efficiently find the shortest path between two points in a graph. We can thus approximate the shortest path between two points in the terrain by running Dijkstra's algorithm, using these points and the terrain's graph representation as input. To increase the accuracy of our answer, we can add more vertices and edges. We construct a graph following the approach of Sack and Lanthier. Let G be the graph from the terrain representation, and let G' be augmented graph that we construct. We add all vertices in G to G'. Between every two vertices that were connected by edges in G, we add a constant number of vertices, equally spaced. These vertices are known as Steiner points. We add an edge between each Steiner point and its neighbor, which will form a triangulated graph. We also connect each of these Steiner points to each of the other Steiner points that are on the two edges that form the opposite side of the triangle. This process will give us a graph that we can use as input to Dijkstra's algorithm. Dijkstra's algorithm will give us the shortest path on our terrain model between the two points. The result we get becomes more accurate as we add more Steiner points to each edge. This accuracy, however, comes at the cost of a larger graph size. An analysis of the size increase is given later in this paper. The process just described gives us a graph in which the edges have weights corresponding to distance. For this constructed graph we use linear interpolation to determine the 3D coordinates of each Steiner point and use the Euclidian distance between the two Steiner points in 3D space as the weight for the edge between the two points. The images linked below show results from an implementation of finding shortest paths using this algorithm: * ShortestPathSampleImages ---+++ Fast Marching One method for computing shortest paths on a surface from raster data is to use fast marching methods. Fast marching methods track the behavior over time of a front that starts from a source point and expands outward at a positive speed. To illustrate, we will first describe the front's behavior on perfectly flat terrain and show how this behavior helps to solve our shortest paths problem. The position of the front at time 0 is our shortest paths source point. We give the front a speed of 1 meter per second. The front expands outward in a circle, and every point in the plane that the front reaches at a time of 1 second is 1 meter from the source. All plane points reached at time 2 are 2 meters away, and so on. The length of the shortest path to any point in the plane is given by the arival time. On a raster surface we will know the elevation of the terrain at regularly spaced sample points. We will call these elevation points; between these elevation points we have a constant number of grid points. All elevation points are grid points, but not all grid points are elevation points. Since the front always expands, it passes each grid point exactly once. Each grid point has a value indicating the arrival time of the front. We can calculate the arrival time at points on a grid using fast marching methods as follows. Kimmel and Sethian tag points as _alive_, _near_, or _far_ [1]. Alive points are points whose arrival times have been determined, while far points are points whose arrival times have not been determined. Near points are points that have an arrival time estimated but do not have the arrival time precisely determined. The algorithm can be summarized as follows: * Find the near point with the smallest arrival time and mark it _alive_. * Mark all neighboring points that are _far_ as _near_. * Update the arrival times of these near points based on the arrival times at neighboring alive points * Repeat if there are still points that are not alive We are able to mark the near point with the smallest arrival time as alive because arrival times are calculated in increasing order and are only calculated from alive points. Initially, the only points that are alive are the points that correspond to the initial position of the front. These points have an arrival time of 0. In our shortest paths problem, the front starts at a single point, so only one point is alive. We update far points adjacent to the alive points and mark them as near. We update the arrival times at near points and repeat until all points are alive. We next describe the procedure for updating the arrival times of points based on arrival times at their alive neighbors. ---++ Updating a grid point's arrival time based on the arrival times of its neighbors The following figure labels the grid points surrounding the grid point to be updated. The arrival time at each grid point is indicated by the T values in the figure. In our examples, at least one of the points a, b, c, or d will be alive, and the point in the center will be near. <img src="%ATTACHURLPATH%/ks_fig2.png" alt="[1]" width="439" height="207" /> [1] _The following needs repair._ Note that _D=R*T_, which you want to use as _T/D = 1/R_, so you actually use reciprocals of speeds... -- Main.JackSnoeyink - 23 May 2005 Suppose that at each grid point _i,j_ we know the front's reciprocal speed \[ F_i,j \] but not the direction; we describe how we find this speed below. We can use the front's speed and the arrival time at one or two adjacent grid points to calculate the arrival time T at grid point _i,j_. We assume that grid points are equally spaced at distance _h_. Decreasing the value of _h_ gives us a more accurate estimates for distance at the cost of more memory. Let p be the grid point whose arrival time we are currently updating. The method we use to update p's arrival time depends how many of p's neighbors are _alive_. If there is only one such neighbor a, let the front's arrival time at a be \[ T_a \]. We assume that the front moves along the horizontal or vertical line between p and a at speed \[ 1/F_i,j \]. The distance between these two grid points is the grid spacing _h_, so the arrival time _T_ at _p_ is _T = Ta + h * Fi,j_: (square both sides) * 1D equation: <br /> <img src="%ATTACHURLPATH%/ks_quad_1d.png" alt="1D equation" width="144" height="30" /> [1] If two of p's neighbors are _alive_ and these two neighbors are not collinear with p, we follow this procedure instead. Let _Ta_ and _Tc_ denote the respective arrival times at adjacent points _a_ and _c_. We use these arrival times in an equation whose unknown is _T_, the arrival time at point _p_. In our problem, we know the speed of the front at point _i,j_. We know the arrival time of the front at grid point _Ta_, which is at a horizontal distance h, and the arrival time at _Tc_, at a vertical distance of h. The front will go a horizontal distance _h_ in time _T-Ta_ and a vertical distance _h_ in time _T-Tc_, giving us a reciprocal velocity vector _( (T-Ta)/h, (T-Tc)/h)_. The length of this vector is the reciprocal of the speed of our front at point i,j. We square both sides to get the quadratic equation below. <img src="%ATTACHURLPATH%/ks_quad.PNG" alt="[1]" /> [1] We can ignore solutions to the quadratic equation that are before the arrival time at _Ta_ or _Tc_, since we know that the front reaches _a_ and _c_ before it reaches _p_. If arrival times at more than two grid points adjacent to the update point p are known, we follow the procedure above for every possible combination of two points where the chosen two points and p are not collinear. We set the arrival time at p to the smallest of these times. This updates the estimated arrival time for near point p. To use this method on a terrain that is not flat we can adjust our speed function in certain places. On flat terrain the speed will still be 1, but when the terrain is not flat the speed can be slower to account for the change in elevation. One easy way to do this is to assign a cost to movement between grid points for which we have the elevation, based on the difference in elevation of these two grid points. We can adjust the speed of the front at intermediate points between these two grid points. The reciprocal of the velocity vector is _((T-Ta)/ sqrt(h^2 + (z(x,y)-z(x-h,y))^2), (T-Tc)/ sqrt(h^2 + (z(x,y)-z(x,y-h))^2))_. That just changes some of the constants in the quadratic. ---++ Evaluating accuracy We use DEM (digital elevation model) sample data taken from an area near Hillsborough, North Carolina. Each DEM is a 500 cell by 500 cell grid, and each cell represents an area that is 20 feet by 20 feet. The height values for cells in the grid are given in feet. We subsample the grid to make the computation faster. We choose two points in the grid, one for the source and one for the destination, and run both algorithms using the same input points. This gives us two shortest paths, one for each algorithm. Each algorithm has a parameter that determines the tradeoff between data size and accuracy. The parameter to the Sack and Lanthier algorithm [2] will be the number of steiner points inserted on each triangle edge, while the parameter to the Kimmel and Sethian [1] algorithm will be the spacing of the grid constructed between known points. ---++ Sack/Lanthier We use a conversion filter to change raster data to the x,y,z points expected by the shortest paths algorithm. We connect the points according to their Delaunay triangulation. We add Steiner points along triangle edges as determined by the input parameter, and we use Dijkstra's algorithm to find the length of the shortest path between the two points. ---++ Kimmel/Sethian The fast marching method works on rasters, so we use the DEM raster data from the source described above. We set the grid spacing as determined by the input parameter. We start the front at the source point, and we use the arrival time at the destination point to find the length of the shortest path. ---++ Evaluation We compare shortest paths from both algorithms as a function of the accuracy parameter (number of Steiner points or grid spacing). The data will be evaluated to determine how changing these parameters affects accuracy and problem size. Specifically, accuracy as a function of input size will be evaluated by using the number of points considered by each algorithm as the input size. For the Sack and Lanthier algorithm, the number of points n, where v and e are the number of vertices and edges in the original graph and s is the number of steiner points added to each edge, is: * n = v + (s * e) For fast marching methods, let r be the number of rows and c be the number of columns in the original grid. The number of vertices v is given by v = r * c. Let h be the number of points added between two grid points. We then have * n = (r + (h * (r - 1)) ) * (c + (h * (c - 1)) ) This works out to * O(n) = O ( h^2 * v + 2 h * v + v ) = O(h^2 * v) Since our implementation makes it difficult to compare input sizes, we use time as our common metric. We give results below of a run on two points. The same two points are used for both algorithms. ---++ Results The following results show how an implementation of Sack and Lanthier's shortest paths algorithm finds a shorter path as the number of Steiner points increases. As the number of points increases, however, so does the number of vertices and edges. |Steiner Points|Distance|Time (ms)| |1|5083.268049|671| |2|4916.525693|942| |3|4879.713764|1692| |4|4875.725698|2594| |5|4869.545668|3856| |6|4867.12319|18497| |7|4864.413176|11037| |8|4863.200293|11146| |9|4861.9047|13510| |10|4861.234855|18137| |11|4860.58039|69653| |12|4860.345447|62342| We can compare this to fast marching method results: |Cell spacing|Distance|Time (ms)| |2|4799.290676|291| |3|4774.699678|971| |4|4761.025498|2694| |5|4752.200848|6189| |6|4745.986264|14512| |7|4741.349517|56613| |8|4737.744565|125706| |9|4734.853847|230010| |10|4732.479382|369536| |11|4730.490979|564945| |12|4728.799324|816216| |13|4727.340999|1141056| |14|4726.069669|1557752| We can see how the difference between the returned shortest path length and the lower bound (Euclidian distance) decreases as our amount of computation time increases. * Comparison of shortest path algorithms: <br /> <img src="%ATTACHURLPATH%/sp_comparison_plot_sm.png" alt="Comparison of shortest path algorithms" width="600" height="450" /> * More comparison: <br /> <img src="%ATTACHURLPATH%/sp_comp_1a.png" alt="More comparison" width="528" height="411" /> * More comparison: <br /> <img src="%ATTACHURLPATH%/sp_comp_2a.png" alt="More comparison" width="533" height="414" /> Of our two implementations, our Sack and Lanthier implementation is limited more by memory than our fast marching implementation. We ran both algorithms with the same amount of Java heap space (a maximum of 384MB), but we ran out of memory after adding 12 Steiner points to each edge. The fast marching algorithm showed a quadratic increase in time, eventually making it impractical to continue. * ShortestPathResults -- more results (raw format) ---++ Implementation bugs Our fast marching algorithm occasionally returns inconsistent results, mainly around the edges of the terrain. This may have to do with our interpolation methods; the algorithm may try to interpolate based on values that lie beyond the edge of the terrain. Our results show only points between shortest paths whose returned values were consistent (the shortest path length decreased as the input size increased). ---++ Alternate approaches The fast marching method could be implemented by performing an 8-way update instead of a 4-way update. We take the arrival time at two adjacent grid points; in this example we have arrival times for Ta and Te. * Figure derived from figure in [1] showing 8-way update: <br /> <img src="%ATTACHURLPATH%/ks_fig2a.png" alt="Figure derived from figure in [1] showing 8-way update" width="439" height="207" /> We know the front covers a horizontal distance of h in time T-Ta, and it covers a distance of h * sqrt(2) that is half horizontal and half vertical in time T-Te. Let k = (T-Tb)/(h*sqrt(2)). We let vector F = [ T-Ta/h + (1/2)k, k ]. The length of this vector gives us the distance the front covers in one time unit, so |F|/1 = |F| is our speed at point i,j. ---++ References * [1] Kimmel, R. and Sethian, J. A., Computing geodesic paths on manifolds, _PNAS_ 95:15, pp. 8431-8435, 1998. * [2] Lanthier, M., Maheshwari, A., and Sack, J-R., Shortest Anisotropic Paths on Terrains, _ICAL '99: Proceedings of the 26th International Colloquium on Automata, Languages and Programming_, Springer-Verlag: London, 1999, pp. 524-533. -- Main.HenryMcEuen - 27 Apr 2005
to top
End of topic
Skip to action links
|
Back to top
Edit
|
Attach image or document
|
Printable version
|
Raw text
|
More topic actions
Revisions: | r1.1
|
Total page history
|
Backlinks
You are here:
TModeling
>
ShortestPathComparison
to top
Copyright © 1999-2024 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki?
Send feedback