Ridge-Valley Path Planning for 3D Terrains - Semantic Scholar
Proceedings the 2006 IEEE International Conference on Robotics and Automation D. Page, A.ofKoschan, M. Abidi, and J. Overholt, "Ridge-Valley Path Planning for 3D Terrains,” in Proc. IEEE International Conference on 221 Orlando, Florida - May 2006 Robotics and Automation ICRA06, Orlando, FL, pp. 119-124, May 2006.
Ridge-Valley Path Planning for 3D Terrains D. L. Page, A. F. Koschan, M. A. Abidi
J. L. Overholt
Imaging, Robotics, and Intelligent Systems The University of Tennessee Ferris Hall 336 Knoxville, Tennessee 37996–2100 Email: {dpage,akoschan,abidi}@utk.edu
U.S. Army TARDEC – Robotics Mobility Laboratory 6501 E. 11 Mile Rd. AMSRD-TAR-R/263 Warren, Michigan 48397–5000 Email: [email protected]
Abstract— This paper presents a tactical path planning algorithm for following ridges or valleys across a 3D terrain. The intent is to generate a path that enables an unmanned vehicle to surveil with maximum observability by traversing the ridges of a terrain or to operate with maximum covertness by navigating the valleys. The input to the algorithm is a 3D triangle mesh model for the terrain of interest. This mesh may be non-uniform and non-regular. Thus, the algorithm leverages research from computer graphics and computer vision to identify ridge-valley features on the terrain. These features serve as “obstacles” for an artificial potential field algorithm. The valleys are obstacles for a surveillance path, or the ridges are obstacles for a covert path. We incorporate geodesic—rather than Euclidean—distances into the potential field formulation to extend path planning to 3D surfaces. We present the theory of our proposed algorithm and provide experimental results.
I. I NTRODUCTION Recent research in 3D computer vision and computer graphics has focused on estimating ridge-valley features over triangle mesh models [1]–[5]. The basic principle of these papers is first to robustly estimate surface curvature and then to identify the ridges and valleys as lines of curvature extrema. The examples in Figs. 1 and 7 illustrate this concept. The terrain maps in these example are triangle meshes consisting of over 40, 000 triangles. We propose to use the ridge features (and conversely the valley features) as obstacles for an artificial potential field algorithm [6] to generate a tactical path through the valleys (or across the ridges) of the terrain. The aerial view in Fig. 7(b) shows such a path from the green sphere to the blue sphere with a closer view of the green starting location in Fig. 1. The block diagram in Fig. 2 outlines the procedures for generating the path. The input to the algorithm is a triangle mesh model of a 3D terrain. One source for such models is digital elevation maps (DEMs). Although most path planning algorithms are formulated with a regular grid structure, these terrain meshes are not necessarily regular nor uniform [7]. The contributions of the paper are as follows: •
•
Formulation of a ridge-valley algorithm for tactical path planning over 3D terrains with non-uniform, non-regular triangulations, Incorporation of geodesic distance as support for potential field construction, and
0-7803-9505-0/06/$20.00 ©2006 IEEE
(a)
(b)
Fig. 1. The proposed algorithm follows the valleys of a terrain to maximize covertness. (a) A close-in view of the terrain highlights the ridges and valleys along with the path of interest. (b) The triangulated surface requires computational methods from computer vision and computer graphics to govern the algorithm.
Integration of 3D geometric processing methods into path planning research. We explore these contributions further relative to other research in Section II. We follow this discussion with Section III, which formulates our proposed algorithm. Then, Section IV presents experimental results. Finally, we conclude in Section V with a summary discussion. •
II. R ELATED W ORK The work that has inspired the development of our algorithm is a series of papers by Lee et al. [8], [9]. The application domain for these papers is not directly related to robotics but is more closely associated with computer graphics. The particular application that Lee et al. target is mesh “scissoring.” This terminology derives from their objective to “cut” a mesh into parts. They apply their scissoring algorithm to 3D boundary meshes that define geometric objects such as horses, bunnies, chairs, and other shapes of interest to the computer graphics community. The basic steps of their algorithm are to select a feature—in their case part boundaries, but for us ridges and valleys—and then threshold to create an initial set of feature lines. This initial set is disjoint and does not form a contiguous line. Thus, Lee et al. define three cost functions over the mesh and then use a weighted average of these three costs to form a single cost. Using a shortest path search on this cost, the scissoring algorithm marches from one end point to
119
Fig. 2. Block diagram of the proposed path planning algorithm for 3D terrains. The input to the algorithm is a DEM triangle mesh, which is possibly non-uniform and non-regular. The output is a path following either the valleys or ridges of the terrain depending on the tactical intent. The attractive field is based solely on the location of the destination—and not the shape of the mesh. (a)
the next of the disjoint feature lines in a tip-to-tail fashion. The result is a continuous unbroken feature line across the surface of the mesh. Although Lee et al. do not directly draw a connection to robotics and path planning, their work has noted similarities. We propose that reformulating Lee et al. to the general framework of path planning avoids the heuristic definitions for their cost function. Another algorithm that has inspired our approach is the fast marching methods of Kimmel and Sethian [10]. They have proposed an algorithm for computing geodesic distances over 3D triangle meshes whereby they “march” from a source point on the mesh outwards. They control the march using a priority queue to insure that the algorithm moves to the next closest vertex in sequence. Given a destination point of interest, their algorithm backtracks along the gradient (of the distance field from the marching process) towards the original source and thus establishes the geodesic path from the destination to the source. As with Lee et al., this method also has a noted similarity to robotics and path planning, but again the connection has not been made in the literature. Although Lee et al. [8], [9] and Kimmel and Sethian [10] have inspired the development of our algorithm, the other area of related work is in the ridge-valley features, themselves. As noted previously, the computer vision and computer graphics communities offer several papers [1]–[5] that address this issue. The one most applicable to our approach is the work of Ohtake et al. [1]. They define a series of equations and an algorithm for extracting ridge-valley feature lines on triangle meshes. These feature lines serve as the starting point for our work. III. T HEORY Following [6], [11]–[13], the basic idea of our algorithm is to create two artificial potential fields, Uatt and Urep over the surface of the triangle mesh and for these fields to govern the path planning. The first field Uatt is an attractive field that pulls the path away from the source location towards the destination. The second field Urep is a repulsive field that pushes the path away from obstacles. In our case, the obstacles are either ridges or valleys of the terrain depending on our tactical intent. The images in Fig. 3 illustrate both of these fields for the path in Fig. 7(b). A “hot” false coloring scheme represents the potential field values such that yellow and red bands are
(b) Fig. 3. Two artificial potential fields govern the path planning algorithm. (a) The first one attracts the path towards the destination. (b) The second one repels the path away from obstacles, in this case ridges. Geodesic distances underly each of the fields rather than Euclidean distances.
closer to the origin and darker black bands are farther away. The alternating white bands further aid visualization. A. Artificial Potential Fields Again, following the path planning literature, if we combine these potential fields (U = Uatt + Urep ), then we can construct a “force” to drive the path from source to destination. The field exerts such a force F as follows: F = −∇U = Fatt + Frep ,
(1)
where ∇ denotes the vector gradient, and Fatt and Frep are the resulting attractive and repulsive forces. Our next task is to define the structure of these artificial potential fields that govern (1). In the discussion below, we use ρgoal (q) to denote a distance function—not necessarily Euclidean—from some point q on the surface of mesh M to some target point qgoal also on M . ρgoal (q) = ρ(q, qgoal ) . Similarly, the distance ρ(q) is a distance function from an “obstacle” on M to some point q also on M . ρ(q) = ρ(q, p) where p is the closest point of the obstacle to q.
120
in Fig. 3(a) illustrates the application of this method to the terrain model. The color-coding scheme alternates bands of white stripes and hot stripes to depict the growing geodesic distance away from the source point—a yellow circle from which the stripes seem to radiate. Similarly, Fig. 3(b) shows geodesic distances from ridge lines of the terrain. C. Ridge-Valley Feature Fig. 4. This sketch demonstrates the difference between a Euclidean distance and a geodesic one. The Euclidean spans the 3D space between the two lobes of the 3D object, but the geodesic distance must remain within the manifold surface that bounds the 3D volume and travel around the lobes.
With these definitions, we select the following function for the attractive potential ⎧ 1 2 ⎨ 2 ζρgoal (q) if ρgoal (q) ≤ d Uatt (q) = , (2) ⎩ dζρgoal (q) if ρgoal (q) > d where ζ is a scaling factor and d is a transition distance. The above function is a hybrid of a parabolic well and conic well to provide stability throughout, but the function is entirely dependent on ρgoal . The repulsive potential—as the name implies—repels the path away from obstacles such that the path never “hits” an obstacle, but unlike the attractive potential, this repulsive potential should have more of a local effect such that obstacles far away have little influence on the path. If we define some distance ρo as the region of influence for an obstacle, then we have the following potential ⎧ 2 1 ⎪ ⎨ 12 η ρ(q) − ρ1o if ρ(q) ≤ ρo , (3) Urep (q) = ⎪ ⎩ 0 if ρ(q) > ρo where η is a scaling constant similar to ζ, and ρo denotes a region of influence. As with d above, the algorithm is not sensitive to ρo , but η has a similar importance as ζ. B. Geodesic Distance We have mentioned distance functions ρgoal (q) and ρ(q) above but without precise definition. In traditional path planning algorithms, these distance functions are the Euclidean distances where, given points q and p, ρ(q, p) = p−q. This definition, however, is not suitable for the distances over our surface mesh M . We instead need a geodesic definition [14] that operates over a discrete surface such as M . Figure 4 shows a visual representation of the difference. The peanut-shaped surface represents a 3D object, and the two curves demonstrate the Euclidean and geodesic distance. The practical application of geodesic distance over a 3D surface to the path planning problem is a proposed contribution of this paper. The fast and efficient computation of geodesics over a discrete surface is not straightforward with a variety of solutions appearing in the literature [10], [15], [16]. For this paper, we use a variation of [16] as an approximate method. The example
Using the geodesic distances above, the next question that we need to consider is the “goal” and “obstacles” that generate the respective Uatt and Urep potentials. The goal, much like the source, is a user designated point on the mesh. We assume that both the goal (i.e. destination) and source locations are given. The field Uatt is the geodesic distance from the destination as in Fig. 3(a). For the obstacles, we must find the ridges (or valleys depending on the tactic) of the terrain mesh. (As an aside, ridges and valleys are the inverse geometrically of one another. For concise presentation, we only consider ridges as obstacles for the remainder of this paper.) If we follow the discussion in [1] for ridges, suppose that our mesh M approximates, at least locally, some smooth surface. Then, at a point on this surface, we can denote the maximum and minimum principal curvatures as kmax and kmin (kmax ≥ kmin ). These scalar values have corresponding principal directions, which we specify as unit vectors Tmax and Tmin . We can further define the derivatives of the principal curvatures along their associated directions as ∂kmax cmax = (4) ∂Tmax where we define cmin in a similar manner. With these derivatives, we can find the extrema of the principal curvatures along their curvature directions with the zero crossings of cmax and cmin . The following conditions characterize ridges: cmax = 0,
∂cmax ∂Tmax
> 0,
kmax > |kmin |.
(5)
Using the above equations and the algorithm in [1], we find the ridge lines as shown in Fig. 7(a). Then using the geodesic distance algorithm in [16], we compute Urep from the ridge lines. The example in Fig. 3(b) shows the results. D. Curvature Estimation To implement (5), we need to estimate principal curvatures and directions along with their associated derivatives. The difficulty is that our mesh M is a singular surface, but several papers [17]–[20] address this issue with the assumption that M approximates a piecewise smooth surface. We follow [20], which is a local surface fitting method. This local fit includes the geometry of the vertices in a one-ring neighborhood and their associated surface normals. The coordinate system of the local parameterization has its origin at the vertex v, and the zaxis aligns with the normal N at v. The x- and y-axis lie in the tangent plane. As with [20], we fit the height function f (x, y) =
121
1 2 1 Ax +Bxy+ Cy 2 +Dx3 +Ex2 y+F xy 2 +Gy 3 , 2 2
which leads to the Weingarten matrix A B −1 W = I II = C D where I and II are the first and second fundamental forms. With eigenanalysis via a Givens rotation, we find the rotation angle θ and thus k1 k2
= A − B tan θ = C + B tan θ
(6)
where we assign kmin and kmax appropriately. We find the principal directions as T1 T2
= X cos θ − Y sin θ = X sin θ + Y cos θ
(a)
(7)
where X and Y are unit-length coordinate vectors in the tangent plane and the assignment to Tmin and Tmax follows (6). For the derivatives of curvature along the principal directions, we follow [21] such that the rank-three derivative tensor is C = ( Dx II Dy II ) where Dx and Dy are derivative operators in the X and Y directions. If we need the derivative of curvature along a specific direction, Tmax for example, we have the following cmax = C(Tmax , Tmax , Tmax ) where the right hand operation denotes multiplying the tensor C by the vector Tmax three times. The resulting scalar cmax is the value that we seek in (4). We show example results for curvature estimation in Fig. 5 for the terrain.
(b) Fig. 5. Curvature estimation for the terrain. (a) At each vertex of the model, we plot the principal directions (red for Tmin , blue for Tmax ). We scale the length of these vectors with the associated principal curvatures (k min , kmax ). (b) A closer view of the ridge of the terrain better illustrates the curvatures.
E. Update Iteration At this point, we have the potential fields Uatt and Urep by estimating surface curvatures, finding the ridges, and then computing geodesic distances. We next start at the source location and walk to the destination under the influence of Fatt and Frep . Several techniques are possible to find these forces. We can use a gradient descent approach, also known as steepest descent or depth-first, as [6] proposed, or we can use other techniques such as best-first search or random search methods. For the remainder of this section, we assume a gradient descent approach. We make this assumption for consistent notation and clear presentation, but the other methods are equally suitable. The following equation defines the subsequent step update: q(k) = q(k−1) + ω
F(k−1) , F(k−1)
(8)
where the superscript (k) denotes the current iteration and ω is the step size. The forces at each iteration are from (1). We begin with (k = 0) and q(0) = qsrc where qsrc is the starting location. We continue until q(k) = qgoal . The step size ω should be small enough so that updates do not jump over an obstacle or overshoot the goal, but ω should not be so small that the path is easily trapped in local minimum. We have found that a practical approach is to tie ω to the average edge length ρe of triangles in M . The
value ρe serves as a uniform scaling parameter. A reasonable value is ω = 0.1ρe . The step size can also vary with the iteration such that ω (k) has a schedule, but we do not utilize this approach. We use a fixed ω throughout the iterations. To implement (8), we must compute the gradients of the potential fields where we recall (1). Locally, if we consider a triangle as in Fig. 6, we estimate the potential at point q within the given triangle by fitting a second-order polynomial. 1 1 Ag x2 + Bg xy + Cg y 2 + Dg x + Eg y + Fg . 2 2 As in [10], the three neighboring triangles support this fitting process such that the potentials at each triangle vertex are the z = g(x, y) values of the function and the x and y values are defined by the local parameterization. The potential gradient thus becomes g(x, y) =
F = (Ag x + Bg y + Dg )X + (Bg x + Cg y + Eg )Y.
(9)
This force controls the iterations and updates of q within a triangle. IV. E XPERIMENTAL R ESULTS With the theory above, we present experimental results in this section. The first results appear in Fig. 7. We derive
122
(a) Fig. 6. The update iterations occur within a triangle. The local fitting procedure defines a local coordinate frame as shown.
(a)
(b) Fig. 7. The 3D terrain model serves as input to our tactical path planning algorithm. (a) The red lines denote the ridges of the terrain that serve as obstacles, in this case, for path planning. (b) The single red line indicates the generated path from the source location (green sphere) at the top of image to the destination (blue sphere) in the lower left.
the terrain in this figure from a DEM model using a simple transformation to form a triangle mesh. (We have downloaded the DEM from the U.S. Geological Survey website at www.usgs.gov.) We apply some mesh smoothing directly to the mesh following a standard Gaussian algorithm. This smoothing aids the curvature estimation algorithms. We first identify the ridge features as in Fig. 7(a). Then we generate the potential fields as in Fig. 3. Applying the iterative gradient descent in Eq. (8), we create the valleyfollowing path in Fig. 7(b). This path begins at the green sphere in the upper right of the figure and continues along the valleys until reaching the blue sphere at the bottom left of the figure. The path faithfully avoids the ridges of the terrain at least until nearing the destination. As the path nears the goal, the transition of the attractive field from a conic to a parabolic well forces the path towards the destination despite
(b)
Fig. 8. The path (red line) traverses over the triangle structure (black lines) of the surface. This triangulation is non-uniform and non-regular. (a) The vector direction from the starting location to the destination defines the initial triangle to begin iterations. (b) The path is free to move across the triangles such that we do not restrict the path to just vertices and edges.
the ridges. The result is that the path begins to prefer a more direct route to the goal as a tradeoff to avoiding the ridges. For this result, we select the parameters as follows. We set the step iteration ω = 0.1ρe as noted previously. The scaling constants are ζ = 0.65 and η = 0.35 to give more importance to the attractive field. The choice for the transition distance is arbitrarily set at d = 10ρe and the region of influence at ρo = 4ρe . The user selects the starting and ending points for the path, but we use the direction vector between these points to establish the initial triangle to begin iterating the path. We project this direction vector into the one-ring neighborhood of the starting vertex and use the local geometry of Fig. 6. To terminate the iterations, we check if the path has entered the one-ring neighborhood of the destination vertex. The close-in views in Fig. 8 illustrate the path (red line) over the triangle structure (black lines). These non-uniform, non-regular triangulations, which these figures emphasize, are different from the regular grid structures often found in 2D path planning. The gap example in Fig. 9(a) is another experimental result. This terrain is a synthetic model that consists of a flat region with a single ridge running across the central section. This ridge has a gap in the middle that bisects the ridge. The triangle mesh consists of 4, 800 triangles. If we select a starting location on one side of the gap and a destination on the other side, the proposed algorithm generates a path that navigates through the gap. This result appears in Fig. 9(b) along with the hot color-coding scheme for the repulsive field. The path is overlaid on the color coding such that the source is in the lower left of the figure and the destination is in the upper right. The parameters are set as ω = 0.1ρe , d = 10ρe , and ρo = 4ρe as in the previous result, but the scaling factors are set slightly different at ζ = 0.8 and η = 0.2. The choice for ζ and η requires careful selection from the user. We have constrained the choice such that ζ +η =1.
(10)
Using the gap example, we see the effects of these parameters in Fig. 10. Suppose we set the attractive field to a small number such that ζ → 0. Then, the attractive field has negligible influence and the repulsive field dominates. The path moves away from the ridge and as a result the gap too as in Fig. 10(a). On the other hand, if we let η → 0, then the opposite
123
avoidance or other maneuvers lead to significant deviations from the original plan. Such an application is the motivation for the work presented in this paper. ACKNOWLEDGMENTS
(a)
(b)
Fig. 9. Gap example. (a) This model consists of a ridge dividing a plain into halves. The ridge has a small gap in the middle that joins these halves together. (b) The proposed algorithm plans a path through the gap in this synthetic terrain model. The repulsive field drives the path away from the the ridges. The hot color-coding scheme illustrates the repulsive potential. The alternating white bands from previous figures are not shown with this visualization since the triangle density is low for this model and such bands would only confuse the visualization.
(a)
(b)
(c)
Fig. 10. With the ridge-gap example, we demonstrate the extremes of the ζ and η parameters. (a) With ζ → 0, the repulsive force dominates the path planning. (b) With η → 0, the attractive force dominates. (c) The desired path occurs with ζ = 0.8 and η = 0.2.
effect occurs. The attractive field dominates and pulls the path straight towards the destination as in Fig. 10(b). The paths does not attempt to avoid the ridges. The desired path appears in Fig. 10(c) where we have set the weights appropriately. V. S UMMARY We have introduced a novel path planning algorithm for unmanned vehicles that allows a tactical plan for following valleys or ridges of a terrain. We have proposed using 3D computational geometry techniques from computer vision and computer graphics to identify the ridges and valleys on nonuniform and non-regular mesh models of terrains. We use these features as obstacles in an potential field path planning algorithm where we have incorporated geodesic distances into this formulation to account for the 3D topology of the terrains. The targeted application for this work is path planning with DEM terrain models. These DEMs serve as input to the algorithm and a tactical path plan is the output. These path plans are suitable for both unmanned ground and unmanned aerial vehicles. One tactic might be for a ground vehicle that needs to covertly travel from a drop zone to some target location. The vehicle would begin by localizing its initial position relative to the DEM. Then, with the coordinates of the target, the proposed algorithm would generate a path that follows the valleys of the terrain. After completing the path plan, the vehicle sets off towards the target. As the vehicle continues toward the target, it constantly localizes position and employs obstacle avoidance. When necessary, the vehicle uses the proposed algorithm again to recompute a path if obstacle
This work is supported by the University Research Program in Robotics under grant DOE-DE-FG52-2004NA25589 and by the DOD/RDECOM/NAC/ARC Program under grant W56HZV-04-2-0001. R EFERENCES [1] Y. Ohtake, A. Belyaev, and H.-P. Seidel, “Ridge-valley lines on meshes via implicit surface fitting,” in Computer Graphics Proceedings (SIGGRAPH 2004), vol. 23, no. 3, 2004, pp. 609–612. [2] V. Interrante, H. Fuchs, and S. Pizer, “Enhancing transparent skin surfaces with ridge and valley lines,” in Proceedings Visualization 1995, 1995, pp. 221–228. [3] G. Stylianou and G. Farin, “Crest lines for surface segmentation and flattening,” IEEE Trans. Visual. Comput. Graphics, vol. 10, no. 5, pp. 536–544, 2004. [4] A. Hubeli and M. Gross, “Multiresolution feature extraction for unstructured meshes,” in Proceedings Visualization 2001, 2001, pp. 287–294. [5] K. Watanabe and A. G. Belyaev, “Detection of salient curvature features on polygonal surfaces,” in Computer Graphics Forum (Eurographics 2001), A. Chalmers and T.-M. Rhyne, Eds., vol. 20, no. 3, 2001, pp. 385–392. [6] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” International Journal of Robotics Research, vol. 5, no. 1, pp. 490–496, Jan. 1986. [7] R. J. Campbell and P. J. Flynn, “A survey of free-form object representation and recognition techniques,” Computer Vision and Image Understanding, vol. 81, pp. 166–201, Nov. 2001. [8] Y. Lee, S. Lee, A. Shamir, D. Cohen-Or, and H.-P. Seidel, “Intellligent mesh scissoring using 3D snakes,” in Proceedings of the Conference on Pacific Graphics, Seoul, Korea, Oct. 2004, pp. 279–287. [9] ——, “Mesh scissoring with minima rule and part salience,” Computer Aided Geometric Design, vol. 22, no. 5, pp. 444–465, July 2005. [10] R. Kimmel and J. A. Sethian, “Computing geodesic paths on manifolds,” in Proceedings of the National Acadmeny of Sciences, vol. 95, no. 15, 1998, pp. 8431–8435. [11] J. Borenstein and Y. Koren, “Real-time obstacle avoidance for fast mobile robots,” IEEE Trans. Syst., Man, Cybern. B, vol. 19, no. 5, pp. 1179–1187, 1989. [12] J. C. Latombe, Robot Motion Planning. Boston, MA: Kluwer Academic Publishers, 1991. [13] S. S. Ge and Y. J. Cui, “New potential functions for robot path planning,” IEEE Trans. Robot. Automat., vol. 16, no. 5, pp. 615–620, Oct. 2000. [14] B. O’Neill, Elementrary Differential Geometry, 2nd ed. Orlando, FL: Academic Press, 1997. [15] K. Polthier and M. Schmies, “Straightest geodesics on polyhedral surfaces,” in Mathematical Visualization, K. P. H. C. Hege, Ed. SpringerVerlag, 1998, pp. 391–409. [16] M. Novotni and R. Klein, “Computing geodesic distances on triangular meshes,” in Proceedings of the WSCG 2002, 2002, pp. 341–347. [17] G. Taubin, “Estimating the tensor of curvature of a surface from a polyhedral approximation,” in Proceedings of the Fifth International Conference on Computer Vision, 1995, pp. 902–907. [18] D. Cohen-Steiner and J.-M. Morvan, “Restricted delaunay triangulations and normal cyle,” in Proceedings ACM Symposium on Computational Geometry 2003, 2003, pp. 237–246. [19] M. Meyer, P. Desbrun, M. Schr¨oder, and A. H. Barr, “Discrete differential geometry operators for triangulated 2-manifolds,” in Proceedings Visualization on Mathematics 2002, H.-C. Hege and K. Polthier, Eds., 2002, pp. 35–57. [20] J. Goldfeather and V. Interrante, “A novel cubic-order algorithm for approximating principal direction vectors,” ACM Transactions on Graphics, vol. 23, no. 1, pp. 45–63, Jan. 2004. [21] S. Rusinkiewicz, “Estimating curvatures and their derivatives on triangle meshes,” in Proceedings of the International Symposium on 3D Data Processing, Visualization, and Transmission, Sept. 2004.
124
Ridge-Valley Path Planning for 3D Terrains D. L. Page, A. F. Koschan, M. A. Abidi
J. L. Overholt
Imaging, Robotics, and Intelligent Systems The University of Tennessee Ferris Hall 336 Knoxville, Tennessee 37996–2100 Email: {dpage,akoschan,abidi}@utk.edu
U.S. Army TARDEC – Robotics Mobility Laboratory 6501 E. 11 Mile Rd. AMSRD-TAR-R/263 Warren, Michigan 48397–5000 Email: [email protected]
Abstract— This paper presents a tactical path planning algorithm for following ridges or valleys across a 3D terrain. The intent is to generate a path that enables an unmanned vehicle to surveil with maximum observability by traversing the ridges of a terrain or to operate with maximum covertness by navigating the valleys. The input to the algorithm is a 3D triangle mesh model for the terrain of interest. This mesh may be non-uniform and non-regular. Thus, the algorithm leverages research from computer graphics and computer vision to identify ridge-valley features on the terrain. These features serve as “obstacles” for an artificial potential field algorithm. The valleys are obstacles for a surveillance path, or the ridges are obstacles for a covert path. We incorporate geodesic—rather than Euclidean—distances into the potential field formulation to extend path planning to 3D surfaces. We present the theory of our proposed algorithm and provide experimental results.
I. I NTRODUCTION Recent research in 3D computer vision and computer graphics has focused on estimating ridge-valley features over triangle mesh models [1]–[5]. The basic principle of these papers is first to robustly estimate surface curvature and then to identify the ridges and valleys as lines of curvature extrema. The examples in Figs. 1 and 7 illustrate this concept. The terrain maps in these example are triangle meshes consisting of over 40, 000 triangles. We propose to use the ridge features (and conversely the valley features) as obstacles for an artificial potential field algorithm [6] to generate a tactical path through the valleys (or across the ridges) of the terrain. The aerial view in Fig. 7(b) shows such a path from the green sphere to the blue sphere with a closer view of the green starting location in Fig. 1. The block diagram in Fig. 2 outlines the procedures for generating the path. The input to the algorithm is a triangle mesh model of a 3D terrain. One source for such models is digital elevation maps (DEMs). Although most path planning algorithms are formulated with a regular grid structure, these terrain meshes are not necessarily regular nor uniform [7]. The contributions of the paper are as follows: •
•
Formulation of a ridge-valley algorithm for tactical path planning over 3D terrains with non-uniform, non-regular triangulations, Incorporation of geodesic distance as support for potential field construction, and
0-7803-9505-0/06/$20.00 ©2006 IEEE
(a)
(b)
Fig. 1. The proposed algorithm follows the valleys of a terrain to maximize covertness. (a) A close-in view of the terrain highlights the ridges and valleys along with the path of interest. (b) The triangulated surface requires computational methods from computer vision and computer graphics to govern the algorithm.
Integration of 3D geometric processing methods into path planning research. We explore these contributions further relative to other research in Section II. We follow this discussion with Section III, which formulates our proposed algorithm. Then, Section IV presents experimental results. Finally, we conclude in Section V with a summary discussion. •
II. R ELATED W ORK The work that has inspired the development of our algorithm is a series of papers by Lee et al. [8], [9]. The application domain for these papers is not directly related to robotics but is more closely associated with computer graphics. The particular application that Lee et al. target is mesh “scissoring.” This terminology derives from their objective to “cut” a mesh into parts. They apply their scissoring algorithm to 3D boundary meshes that define geometric objects such as horses, bunnies, chairs, and other shapes of interest to the computer graphics community. The basic steps of their algorithm are to select a feature—in their case part boundaries, but for us ridges and valleys—and then threshold to create an initial set of feature lines. This initial set is disjoint and does not form a contiguous line. Thus, Lee et al. define three cost functions over the mesh and then use a weighted average of these three costs to form a single cost. Using a shortest path search on this cost, the scissoring algorithm marches from one end point to
119
Fig. 2. Block diagram of the proposed path planning algorithm for 3D terrains. The input to the algorithm is a DEM triangle mesh, which is possibly non-uniform and non-regular. The output is a path following either the valleys or ridges of the terrain depending on the tactical intent. The attractive field is based solely on the location of the destination—and not the shape of the mesh. (a)
the next of the disjoint feature lines in a tip-to-tail fashion. The result is a continuous unbroken feature line across the surface of the mesh. Although Lee et al. do not directly draw a connection to robotics and path planning, their work has noted similarities. We propose that reformulating Lee et al. to the general framework of path planning avoids the heuristic definitions for their cost function. Another algorithm that has inspired our approach is the fast marching methods of Kimmel and Sethian [10]. They have proposed an algorithm for computing geodesic distances over 3D triangle meshes whereby they “march” from a source point on the mesh outwards. They control the march using a priority queue to insure that the algorithm moves to the next closest vertex in sequence. Given a destination point of interest, their algorithm backtracks along the gradient (of the distance field from the marching process) towards the original source and thus establishes the geodesic path from the destination to the source. As with Lee et al., this method also has a noted similarity to robotics and path planning, but again the connection has not been made in the literature. Although Lee et al. [8], [9] and Kimmel and Sethian [10] have inspired the development of our algorithm, the other area of related work is in the ridge-valley features, themselves. As noted previously, the computer vision and computer graphics communities offer several papers [1]–[5] that address this issue. The one most applicable to our approach is the work of Ohtake et al. [1]. They define a series of equations and an algorithm for extracting ridge-valley feature lines on triangle meshes. These feature lines serve as the starting point for our work. III. T HEORY Following [6], [11]–[13], the basic idea of our algorithm is to create two artificial potential fields, Uatt and Urep over the surface of the triangle mesh and for these fields to govern the path planning. The first field Uatt is an attractive field that pulls the path away from the source location towards the destination. The second field Urep is a repulsive field that pushes the path away from obstacles. In our case, the obstacles are either ridges or valleys of the terrain depending on our tactical intent. The images in Fig. 3 illustrate both of these fields for the path in Fig. 7(b). A “hot” false coloring scheme represents the potential field values such that yellow and red bands are
(b) Fig. 3. Two artificial potential fields govern the path planning algorithm. (a) The first one attracts the path towards the destination. (b) The second one repels the path away from obstacles, in this case ridges. Geodesic distances underly each of the fields rather than Euclidean distances.
closer to the origin and darker black bands are farther away. The alternating white bands further aid visualization. A. Artificial Potential Fields Again, following the path planning literature, if we combine these potential fields (U = Uatt + Urep ), then we can construct a “force” to drive the path from source to destination. The field exerts such a force F as follows: F = −∇U = Fatt + Frep ,
(1)
where ∇ denotes the vector gradient, and Fatt and Frep are the resulting attractive and repulsive forces. Our next task is to define the structure of these artificial potential fields that govern (1). In the discussion below, we use ρgoal (q) to denote a distance function—not necessarily Euclidean—from some point q on the surface of mesh M to some target point qgoal also on M . ρgoal (q) = ρ(q, qgoal ) . Similarly, the distance ρ(q) is a distance function from an “obstacle” on M to some point q also on M . ρ(q) = ρ(q, p) where p is the closest point of the obstacle to q.
120
in Fig. 3(a) illustrates the application of this method to the terrain model. The color-coding scheme alternates bands of white stripes and hot stripes to depict the growing geodesic distance away from the source point—a yellow circle from which the stripes seem to radiate. Similarly, Fig. 3(b) shows geodesic distances from ridge lines of the terrain. C. Ridge-Valley Feature Fig. 4. This sketch demonstrates the difference between a Euclidean distance and a geodesic one. The Euclidean spans the 3D space between the two lobes of the 3D object, but the geodesic distance must remain within the manifold surface that bounds the 3D volume and travel around the lobes.
With these definitions, we select the following function for the attractive potential ⎧ 1 2 ⎨ 2 ζρgoal (q) if ρgoal (q) ≤ d Uatt (q) = , (2) ⎩ dζρgoal (q) if ρgoal (q) > d where ζ is a scaling factor and d is a transition distance. The above function is a hybrid of a parabolic well and conic well to provide stability throughout, but the function is entirely dependent on ρgoal . The repulsive potential—as the name implies—repels the path away from obstacles such that the path never “hits” an obstacle, but unlike the attractive potential, this repulsive potential should have more of a local effect such that obstacles far away have little influence on the path. If we define some distance ρo as the region of influence for an obstacle, then we have the following potential ⎧ 2 1 ⎪ ⎨ 12 η ρ(q) − ρ1o if ρ(q) ≤ ρo , (3) Urep (q) = ⎪ ⎩ 0 if ρ(q) > ρo where η is a scaling constant similar to ζ, and ρo denotes a region of influence. As with d above, the algorithm is not sensitive to ρo , but η has a similar importance as ζ. B. Geodesic Distance We have mentioned distance functions ρgoal (q) and ρ(q) above but without precise definition. In traditional path planning algorithms, these distance functions are the Euclidean distances where, given points q and p, ρ(q, p) = p−q. This definition, however, is not suitable for the distances over our surface mesh M . We instead need a geodesic definition [14] that operates over a discrete surface such as M . Figure 4 shows a visual representation of the difference. The peanut-shaped surface represents a 3D object, and the two curves demonstrate the Euclidean and geodesic distance. The practical application of geodesic distance over a 3D surface to the path planning problem is a proposed contribution of this paper. The fast and efficient computation of geodesics over a discrete surface is not straightforward with a variety of solutions appearing in the literature [10], [15], [16]. For this paper, we use a variation of [16] as an approximate method. The example
Using the geodesic distances above, the next question that we need to consider is the “goal” and “obstacles” that generate the respective Uatt and Urep potentials. The goal, much like the source, is a user designated point on the mesh. We assume that both the goal (i.e. destination) and source locations are given. The field Uatt is the geodesic distance from the destination as in Fig. 3(a). For the obstacles, we must find the ridges (or valleys depending on the tactic) of the terrain mesh. (As an aside, ridges and valleys are the inverse geometrically of one another. For concise presentation, we only consider ridges as obstacles for the remainder of this paper.) If we follow the discussion in [1] for ridges, suppose that our mesh M approximates, at least locally, some smooth surface. Then, at a point on this surface, we can denote the maximum and minimum principal curvatures as kmax and kmin (kmax ≥ kmin ). These scalar values have corresponding principal directions, which we specify as unit vectors Tmax and Tmin . We can further define the derivatives of the principal curvatures along their associated directions as ∂kmax cmax = (4) ∂Tmax where we define cmin in a similar manner. With these derivatives, we can find the extrema of the principal curvatures along their curvature directions with the zero crossings of cmax and cmin . The following conditions characterize ridges: cmax = 0,
∂cmax ∂Tmax
> 0,
kmax > |kmin |.
(5)
Using the above equations and the algorithm in [1], we find the ridge lines as shown in Fig. 7(a). Then using the geodesic distance algorithm in [16], we compute Urep from the ridge lines. The example in Fig. 3(b) shows the results. D. Curvature Estimation To implement (5), we need to estimate principal curvatures and directions along with their associated derivatives. The difficulty is that our mesh M is a singular surface, but several papers [17]–[20] address this issue with the assumption that M approximates a piecewise smooth surface. We follow [20], which is a local surface fitting method. This local fit includes the geometry of the vertices in a one-ring neighborhood and their associated surface normals. The coordinate system of the local parameterization has its origin at the vertex v, and the zaxis aligns with the normal N at v. The x- and y-axis lie in the tangent plane. As with [20], we fit the height function f (x, y) =
121
1 2 1 Ax +Bxy+ Cy 2 +Dx3 +Ex2 y+F xy 2 +Gy 3 , 2 2
which leads to the Weingarten matrix A B −1 W = I II = C D where I and II are the first and second fundamental forms. With eigenanalysis via a Givens rotation, we find the rotation angle θ and thus k1 k2
= A − B tan θ = C + B tan θ
(6)
where we assign kmin and kmax appropriately. We find the principal directions as T1 T2
= X cos θ − Y sin θ = X sin θ + Y cos θ
(a)
(7)
where X and Y are unit-length coordinate vectors in the tangent plane and the assignment to Tmin and Tmax follows (6). For the derivatives of curvature along the principal directions, we follow [21] such that the rank-three derivative tensor is C = ( Dx II Dy II ) where Dx and Dy are derivative operators in the X and Y directions. If we need the derivative of curvature along a specific direction, Tmax for example, we have the following cmax = C(Tmax , Tmax , Tmax ) where the right hand operation denotes multiplying the tensor C by the vector Tmax three times. The resulting scalar cmax is the value that we seek in (4). We show example results for curvature estimation in Fig. 5 for the terrain.
(b) Fig. 5. Curvature estimation for the terrain. (a) At each vertex of the model, we plot the principal directions (red for Tmin , blue for Tmax ). We scale the length of these vectors with the associated principal curvatures (k min , kmax ). (b) A closer view of the ridge of the terrain better illustrates the curvatures.
E. Update Iteration At this point, we have the potential fields Uatt and Urep by estimating surface curvatures, finding the ridges, and then computing geodesic distances. We next start at the source location and walk to the destination under the influence of Fatt and Frep . Several techniques are possible to find these forces. We can use a gradient descent approach, also known as steepest descent or depth-first, as [6] proposed, or we can use other techniques such as best-first search or random search methods. For the remainder of this section, we assume a gradient descent approach. We make this assumption for consistent notation and clear presentation, but the other methods are equally suitable. The following equation defines the subsequent step update: q(k) = q(k−1) + ω
F(k−1) , F(k−1)
(8)
where the superscript (k) denotes the current iteration and ω is the step size. The forces at each iteration are from (1). We begin with (k = 0) and q(0) = qsrc where qsrc is the starting location. We continue until q(k) = qgoal . The step size ω should be small enough so that updates do not jump over an obstacle or overshoot the goal, but ω should not be so small that the path is easily trapped in local minimum. We have found that a practical approach is to tie ω to the average edge length ρe of triangles in M . The
value ρe serves as a uniform scaling parameter. A reasonable value is ω = 0.1ρe . The step size can also vary with the iteration such that ω (k) has a schedule, but we do not utilize this approach. We use a fixed ω throughout the iterations. To implement (8), we must compute the gradients of the potential fields where we recall (1). Locally, if we consider a triangle as in Fig. 6, we estimate the potential at point q within the given triangle by fitting a second-order polynomial. 1 1 Ag x2 + Bg xy + Cg y 2 + Dg x + Eg y + Fg . 2 2 As in [10], the three neighboring triangles support this fitting process such that the potentials at each triangle vertex are the z = g(x, y) values of the function and the x and y values are defined by the local parameterization. The potential gradient thus becomes g(x, y) =
F = (Ag x + Bg y + Dg )X + (Bg x + Cg y + Eg )Y.
(9)
This force controls the iterations and updates of q within a triangle. IV. E XPERIMENTAL R ESULTS With the theory above, we present experimental results in this section. The first results appear in Fig. 7. We derive
122
(a) Fig. 6. The update iterations occur within a triangle. The local fitting procedure defines a local coordinate frame as shown.
(a)
(b) Fig. 7. The 3D terrain model serves as input to our tactical path planning algorithm. (a) The red lines denote the ridges of the terrain that serve as obstacles, in this case, for path planning. (b) The single red line indicates the generated path from the source location (green sphere) at the top of image to the destination (blue sphere) in the lower left.
the terrain in this figure from a DEM model using a simple transformation to form a triangle mesh. (We have downloaded the DEM from the U.S. Geological Survey website at www.usgs.gov.) We apply some mesh smoothing directly to the mesh following a standard Gaussian algorithm. This smoothing aids the curvature estimation algorithms. We first identify the ridge features as in Fig. 7(a). Then we generate the potential fields as in Fig. 3. Applying the iterative gradient descent in Eq. (8), we create the valleyfollowing path in Fig. 7(b). This path begins at the green sphere in the upper right of the figure and continues along the valleys until reaching the blue sphere at the bottom left of the figure. The path faithfully avoids the ridges of the terrain at least until nearing the destination. As the path nears the goal, the transition of the attractive field from a conic to a parabolic well forces the path towards the destination despite
(b)
Fig. 8. The path (red line) traverses over the triangle structure (black lines) of the surface. This triangulation is non-uniform and non-regular. (a) The vector direction from the starting location to the destination defines the initial triangle to begin iterations. (b) The path is free to move across the triangles such that we do not restrict the path to just vertices and edges.
the ridges. The result is that the path begins to prefer a more direct route to the goal as a tradeoff to avoiding the ridges. For this result, we select the parameters as follows. We set the step iteration ω = 0.1ρe as noted previously. The scaling constants are ζ = 0.65 and η = 0.35 to give more importance to the attractive field. The choice for the transition distance is arbitrarily set at d = 10ρe and the region of influence at ρo = 4ρe . The user selects the starting and ending points for the path, but we use the direction vector between these points to establish the initial triangle to begin iterating the path. We project this direction vector into the one-ring neighborhood of the starting vertex and use the local geometry of Fig. 6. To terminate the iterations, we check if the path has entered the one-ring neighborhood of the destination vertex. The close-in views in Fig. 8 illustrate the path (red line) over the triangle structure (black lines). These non-uniform, non-regular triangulations, which these figures emphasize, are different from the regular grid structures often found in 2D path planning. The gap example in Fig. 9(a) is another experimental result. This terrain is a synthetic model that consists of a flat region with a single ridge running across the central section. This ridge has a gap in the middle that bisects the ridge. The triangle mesh consists of 4, 800 triangles. If we select a starting location on one side of the gap and a destination on the other side, the proposed algorithm generates a path that navigates through the gap. This result appears in Fig. 9(b) along with the hot color-coding scheme for the repulsive field. The path is overlaid on the color coding such that the source is in the lower left of the figure and the destination is in the upper right. The parameters are set as ω = 0.1ρe , d = 10ρe , and ρo = 4ρe as in the previous result, but the scaling factors are set slightly different at ζ = 0.8 and η = 0.2. The choice for ζ and η requires careful selection from the user. We have constrained the choice such that ζ +η =1.
(10)
Using the gap example, we see the effects of these parameters in Fig. 10. Suppose we set the attractive field to a small number such that ζ → 0. Then, the attractive field has negligible influence and the repulsive field dominates. The path moves away from the ridge and as a result the gap too as in Fig. 10(a). On the other hand, if we let η → 0, then the opposite
123
avoidance or other maneuvers lead to significant deviations from the original plan. Such an application is the motivation for the work presented in this paper. ACKNOWLEDGMENTS
(a)
(b)
Fig. 9. Gap example. (a) This model consists of a ridge dividing a plain into halves. The ridge has a small gap in the middle that joins these halves together. (b) The proposed algorithm plans a path through the gap in this synthetic terrain model. The repulsive field drives the path away from the the ridges. The hot color-coding scheme illustrates the repulsive potential. The alternating white bands from previous figures are not shown with this visualization since the triangle density is low for this model and such bands would only confuse the visualization.
(a)
(b)
(c)
Fig. 10. With the ridge-gap example, we demonstrate the extremes of the ζ and η parameters. (a) With ζ → 0, the repulsive force dominates the path planning. (b) With η → 0, the attractive force dominates. (c) The desired path occurs with ζ = 0.8 and η = 0.2.
effect occurs. The attractive field dominates and pulls the path straight towards the destination as in Fig. 10(b). The paths does not attempt to avoid the ridges. The desired path appears in Fig. 10(c) where we have set the weights appropriately. V. S UMMARY We have introduced a novel path planning algorithm for unmanned vehicles that allows a tactical plan for following valleys or ridges of a terrain. We have proposed using 3D computational geometry techniques from computer vision and computer graphics to identify the ridges and valleys on nonuniform and non-regular mesh models of terrains. We use these features as obstacles in an potential field path planning algorithm where we have incorporated geodesic distances into this formulation to account for the 3D topology of the terrains. The targeted application for this work is path planning with DEM terrain models. These DEMs serve as input to the algorithm and a tactical path plan is the output. These path plans are suitable for both unmanned ground and unmanned aerial vehicles. One tactic might be for a ground vehicle that needs to covertly travel from a drop zone to some target location. The vehicle would begin by localizing its initial position relative to the DEM. Then, with the coordinates of the target, the proposed algorithm would generate a path that follows the valleys of the terrain. After completing the path plan, the vehicle sets off towards the target. As the vehicle continues toward the target, it constantly localizes position and employs obstacle avoidance. When necessary, the vehicle uses the proposed algorithm again to recompute a path if obstacle
This work is supported by the University Research Program in Robotics under grant DOE-DE-FG52-2004NA25589 and by the DOD/RDECOM/NAC/ARC Program under grant W56HZV-04-2-0001. R EFERENCES [1] Y. Ohtake, A. Belyaev, and H.-P. Seidel, “Ridge-valley lines on meshes via implicit surface fitting,” in Computer Graphics Proceedings (SIGGRAPH 2004), vol. 23, no. 3, 2004, pp. 609–612. [2] V. Interrante, H. Fuchs, and S. Pizer, “Enhancing transparent skin surfaces with ridge and valley lines,” in Proceedings Visualization 1995, 1995, pp. 221–228. [3] G. Stylianou and G. Farin, “Crest lines for surface segmentation and flattening,” IEEE Trans. Visual. Comput. Graphics, vol. 10, no. 5, pp. 536–544, 2004. [4] A. Hubeli and M. Gross, “Multiresolution feature extraction for unstructured meshes,” in Proceedings Visualization 2001, 2001, pp. 287–294. [5] K. Watanabe and A. G. Belyaev, “Detection of salient curvature features on polygonal surfaces,” in Computer Graphics Forum (Eurographics 2001), A. Chalmers and T.-M. Rhyne, Eds., vol. 20, no. 3, 2001, pp. 385–392. [6] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” International Journal of Robotics Research, vol. 5, no. 1, pp. 490–496, Jan. 1986. [7] R. J. Campbell and P. J. Flynn, “A survey of free-form object representation and recognition techniques,” Computer Vision and Image Understanding, vol. 81, pp. 166–201, Nov. 2001. [8] Y. Lee, S. Lee, A. Shamir, D. Cohen-Or, and H.-P. Seidel, “Intellligent mesh scissoring using 3D snakes,” in Proceedings of the Conference on Pacific Graphics, Seoul, Korea, Oct. 2004, pp. 279–287. [9] ——, “Mesh scissoring with minima rule and part salience,” Computer Aided Geometric Design, vol. 22, no. 5, pp. 444–465, July 2005. [10] R. Kimmel and J. A. Sethian, “Computing geodesic paths on manifolds,” in Proceedings of the National Acadmeny of Sciences, vol. 95, no. 15, 1998, pp. 8431–8435. [11] J. Borenstein and Y. Koren, “Real-time obstacle avoidance for fast mobile robots,” IEEE Trans. Syst., Man, Cybern. B, vol. 19, no. 5, pp. 1179–1187, 1989. [12] J. C. Latombe, Robot Motion Planning. Boston, MA: Kluwer Academic Publishers, 1991. [13] S. S. Ge and Y. J. Cui, “New potential functions for robot path planning,” IEEE Trans. Robot. Automat., vol. 16, no. 5, pp. 615–620, Oct. 2000. [14] B. O’Neill, Elementrary Differential Geometry, 2nd ed. Orlando, FL: Academic Press, 1997. [15] K. Polthier and M. Schmies, “Straightest geodesics on polyhedral surfaces,” in Mathematical Visualization, K. P. H. C. Hege, Ed. SpringerVerlag, 1998, pp. 391–409. [16] M. Novotni and R. Klein, “Computing geodesic distances on triangular meshes,” in Proceedings of the WSCG 2002, 2002, pp. 341–347. [17] G. Taubin, “Estimating the tensor of curvature of a surface from a polyhedral approximation,” in Proceedings of the Fifth International Conference on Computer Vision, 1995, pp. 902–907. [18] D. Cohen-Steiner and J.-M. Morvan, “Restricted delaunay triangulations and normal cyle,” in Proceedings ACM Symposium on Computational Geometry 2003, 2003, pp. 237–246. [19] M. Meyer, P. Desbrun, M. Schr¨oder, and A. H. Barr, “Discrete differential geometry operators for triangulated 2-manifolds,” in Proceedings Visualization on Mathematics 2002, H.-C. Hege and K. Polthier, Eds., 2002, pp. 35–57. [20] J. Goldfeather and V. Interrante, “A novel cubic-order algorithm for approximating principal direction vectors,” ACM Transactions on Graphics, vol. 23, no. 1, pp. 45–63, Jan. 2004. [21] S. Rusinkiewicz, “Estimating curvatures and their derivatives on triangle meshes,” in Proceedings of the International Symposium on 3D Data Processing, Visualization, and Transmission, Sept. 2004.
124
