# Estimation of road centerline curvature from raw GPS data/ Centrines kelio linijos kreives nustatymas remiantis neapdorotais GPS duomenimis/ Cela ass liekuma noteiksana no neapstradatiem GPS datiem/ Tee telgjoone koveruse hindamine tootlemata GPS-andmete alusel.

1. IntroductionPrecise digital road geometry is a fundamental for development of various applications regarding intelligent transport systems, traffic safety and other traffic-related areas (Barai 2003).

The accuracy of available digital maps suffices for some useful applications, such as navigation and route guidance systems. However, there are other applications, where a combination of enhanced digital road maps and precise positioning systems is necessary. Some possible applications for which detailed road geometry is needed are Rollover Warning and Control system, Lane Departure Warning (LDW) and lane-level navigation. Another application areas that can benefit from detailed road geometry data are analysis of road geometry on driving conditions and visibility (Pellegrino 2009; Vorobjovas 2011) and consistency analysis of highway elements and their influence on operating speed and traffic safety (Cafiso et al. 2010; D'Attoma 2010; Dell'Acqua, Russo 2010; Discetti 2010; Perco 2008; Zuriaga et al. 2010; Sliupas 2009).

There are several approaches to obtain the data. One widely used approach is obtaining the required data with specially equipped cars. The other is the use of the 2D aerial images, photointerpretation and appropriate vectorization. Due to the complexity of such images, the procedure is practically impossible to automate. Another option is statistical approach, where a large quantity of possibly noisy data from global positioning systems (GPS) for a fleet of vehicles is combined. The data are obtained from vehicles that commute on their usual business. This approach can be rather effective and less expensive. The same is true for the price of differential (DGPS) devices. This way one can construct maps that have higher accuracy, are easy to maintain, and are cheaper.

Accurate determination of the road centerline is necessary for deployment of advanced applications, such as LDW for lateral control. These accuracy requirements were integrated within the NextMAP project that deals with the economical and technical aspects of digitally defined road maps. The quality of geometry includes adequate topology of the road network, accuracy of the centerline and last but not least, frequency of the map update. The combination of all these elements guarantees suitability for most demanding applications.

A specially equipped car, as for example Photobus (Gontran et al. 2005), combines accurate positioning through GPS/IMU measurements with a vertically oriented CCD or CMOS camera. The system performs synchronization of GPS data with imagery and allows precise determination of the road centerline. The geodetic coordinates must be transformed to a local reference frame. This data allows the acquisition of the horizontal and vertical components of the road geometry.

The purpose of presented research is to develop and verify an efficient method for curve fitting, which can be used in modeling and analysis of the road centerline axis geometry. To model and efficiently analyze the road geometry, cubic splines were chosen as interpolating functions. Cubic splines are piecewise interpolated curves by 3rd-degree polynomials between n adjustment points. Consequently, they take into account all trajectory behavior on each piece of interval and provide continuity conditions, velocity and acceleration on each adjustment point. Due to these characteristics, such curves are wellsuited for approximating road axes. The next task is to segment given digital curve, representing the road axis, into arcs. It is well known that long straight stretches of road can increase the chance of accidents, as they dull the driver's attention or stimulate him or her to speed up. On the other hand, most turns generate a decrease of speed and can sometimes surprise the driver, especially in case of bad weather conditions. The Advanced Driver Assistance Systems (ADAS) function, the Curve Speed Warning (CSW), or Rollover Warning, can inform users if they are traveling too fast to successfully pass an upcoming turn. Since roadway design guidelines link the radius of turn curvature with the max speed for that turn, information on curvature radii together with other data can make CSW reliable and useful.

2. Overview of published methods

Curve-fitting is an important technique in the processing and analyzing digitaly given curves. Known procedures can be grouped into interpolation and approximation, depending on whether the resulting curve passes through all of the data points or not. To fit a digitaly given curve it has to be divided the into segments. Such digital segment is fitted with a appropriate analytic curve, which can be a line segment, a circular arc, or a high order curve.

The simplest approach, polygonal approximation (Perez, Vidal 1994; Pikaz, Averbuch 1996; Ray, Ray 1993, 1994), fits each digital segment with a straight line segment by finding the breakpoints and connecting them with straight line segments.

Better representation can be obtained by using circular arcs for segment fitting (Pei, Horng 1995; Pei, Horng 1996). Many methods are based on a combination of line segments and circular arcs (Horng, Li 2001; Ichoku et al. 1996; Rosin, West 1989). A significant effort has also been devoted to fitting arcs and bi-arcs through points in 2D (Yang, Du 1996). None of the above methods can be easily extended to 3D cases. Three-dimensional curves play a fundamental role in many 3D applications, one of them is also road centerline in engineering praxis. For this purpose, non-linear curves have been used (Farin 1997; Laurent et al. 1991).

Curvature is one of the most important pieces of information about shape contours, but in spite of its importance, no definitive numerical method for curvature estimation has been found. While the curvature of continuous curves can be precisely evaluated by using closed form expressions, the problem of estimating curvature of spatially sampled digital curves is not straightforward. The main problem with such "digital curvature" estimation approaches is that digitaly given curves do not even have a curvature in a strict sense, for they are no more than a set of isolated points. Some pre-processing steps are needed before "curvature" can be estimated. The importance of knowing this element has motivated a series of researches (Coeurjolly et al. 2001; Estrozi et al. 1999; Fairney, Fairney 1994; Imran et al. 2006; Lewiner et al. 2004; Mokhtarian, Mackworth 1992; Worring, Smeulders 1993).

The procedure for road centerline determination should also estimate the corresponding road curvature. Using a spline-fitting approach, the curvature can be calculated through the 1st and 2nd derivative (Schroedl et al. 2004). It is usually assumed that curvature is continuous across a road segment, which is why the degree of the polynomial used for interpolation must be at least three.

3. The proposed method

Here, a brief description of our approach is presented, which is different from procedures used usually.

The raw GPS data, where points are not equidistant, are smoothed by splines, and new equidistant sample data are generated. These data are used in further computations. The next step is segmentation of digitally given curve into circular arcs. For this purpose some points (10 to 20) were selected and search for the best fitting circular arc was performed. In each further step one additional point is added while observing the error, which is defined as the max distance of chosen points to the arc. When this error exceeds the prescribed error (e.g. 1 m), the construction of the arc is stopped and construction of a new arc starts.

The arc construction is performed in two steps. The abovementioned points are projected stereographically onto a sphere and find the plane whose intersection with the sphere best fits the data points. The intersection of the plane with the sphere is a circle, which is then stereographically projected back onto the XY plane. The result is the corresponding circular arc.

All steps of the presented procedure are numerically very stable and there is no numerical problem when estimating a radius of curvature of nearly linear segments of a digital curve. As a consequence, the proposed method is extremely robust. A detailed outline of the proposed procedure is given below.

The GPS coordinates of road centerline are not equidistant; distances between adjacent points vary from few to several 10 m. That is why the digital curve has to be smoothened.

Let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (1)

and let T be the matrix whose rows are coordinates of points [T.sub.1], [T.sub.2], [T.sub.3], [T.sub.4]. Then

B(f) = [[t.sup.3], [t.sup.2], t, 1] MT (2)

is a cubic polynomial mapping, represented with a matrix product of three matrices of dimensions 1x4, 4x4 and 4x2. While the parameter t takes all values from the interval [0, 1], the point B(t) lies on the arc, starting near point [T.sub.2] and ending near point [T.sub.3]. Then it is necessary to repeat the procedure with points [T.sub.2], [T.sub.3], [T.sub.4], [T.sub.5]. Both polynomial arcs are smoothly ([C.sup.2]) joined near point [T.sub.3] (Bartels et al. 1987). The whole polygonal line can be approximated in that way, except the 1st and the last point. The method works on polygonal lines in arbitrary Euclidean space Rn. The example in [R.sup.2] is illustrated in Fig. 1.

[FIGURE 1 OMITTED]

For the purpose of the problem, a polygonal line whose vertices are equidistant has to be generated. Let [S.sup.2] be the unit sphere with center in the coordinate origin and north pole at N(0, 0, 1). The line through the North pole N and an arbitrary point (u, v, 0) on the equatorial plane intersects the sphere in (x, y, z), if for some real t the following is true:

(x, y, z) = t(u, v, 0) + (1 - t)(0,0,1). (3)

Taking into account that the vector (x, y, z) has unit length, it follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (4)

and

t = 2/[u.sup.2] + [v.sup.2] + 1 (5)

Point (u, v, 0) in the equatorial plane maps into point T(x, y, z) on the sphere, whose coordinates are (Fig. 2)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (6)

[FIGURE 2 OMITTED]

The inverse relationship between point T(x, y, z) on the sphere and point (u, v, 0) on the equatorial plane is

u = x/1-z, v = y/1-z (7)

The lines and circles in the equatorial plane map onto circles on the sphere. The circles passing through the north pole correspond to straight lines on the equatorial plane. To see that one can take an arbitrary circle or line lying on the equatorial (u, v) plane, it can be written as:

{(u,v); a([u.sup.2] + [v.sup.2]) + bu + cv + d = } (8)

and express it in coordinates on the sphere

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

The above expression simplifies into the Eq of the plane

bx + cy + (a-d)z + a + d = 0, (10)

intersecting the sphere in a circle (Fig. 3).

The stereographic projection maps the circle resulting from the intersection of the unit sphere with the plane

Ax + By + Cz = D, (11)

which does not intersect the North pole, onto a circle lying in the equatorial plane with radius

[FIGURE 3 OMITTED]

R = [square root of ([A.sup.2] + [B.sup.2] + [C.sup.2] - [D.sup.2]/[absolute value of C-D] (12)

and center at

S( A/D-C, B/D-C) (13)

For the part L of the above generated polyline a circular arc or line approximation has to be find in the following way. First the stereographically project plane vertices of L onto a sphere must be done. The least square method to fit these 3D space points with the plane [PSI], intersecting the unit sphere in the circle [K.sub.1] is used. Its projection back to the equatorial plane gives the desired circular arc or line [K.sub.2].

For given points

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

there exists such a plane Ax + By + Cz = D that the sum of squares of distances of given points from the plane is min.

First the standard least square problem must be solved, so that a linear affine function f(x,y) = [alpha]x + [beta]y + [gamma], is found - which minimizes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (14)

Then an orthonormal basis with its last base vector orthogonal to the plane is found. The procedure is repeated in the new base. The convergence is extremely fast and only a few steps are needed.

Fig. 4. Mathematica Code baza [a_,b_]; * Module [{}, e1 = a/[square root of (a.a);] e2 = b = (b.e1) e1; e2 = e2 /[square root of (e2.e2)]; {e1,e2, e1 x e2} ] ravnina [d0_] : * Module [n{n,dd}, korak : * Module [{}, f[x_,y_] = Fit [data, {1,x,y}, {x,y}]; h = Coefficient [f[x,y], x]; k = Coefficient [f[x,y], y]; Q = baza [{1, 0, h}, {0, 1, k}] // Transpose; P = Inverse [Q]; data = P.#&/@data; n = {1,0,h} x {0, 1, k}; n = R.n [square root of (n.n)]; dd = n.R. {0, 0, f [0, 0]}; R = R.Q; ]; R = IdentityMatrix[3]; data = d0; Do[korak], {10}]; Sign [dd] {n, dd} //Flatten ];

The Mathematica Code of the above procedure is as follows in Fig. 4.

For testing purposes, an artificial digital curve has been constructed with two circular arcs of radius R = 300 m and R = 200 m and a connecting straight line segment (Fig. 5). The curve points are between 20 m and 35 m apart. The coordinates of the testing curve have been randomly uniformly perturbed inside a square of 2 m by 2 m around each point. This perturbation approx corresponds to the errors in position obtained by a GPS device. The procedure described above yielded two circular arcs of R = 299 m and R = 200 m respectively and a circular arc with R = 64 205 m, which corresponds to the straight line segment connecting both arcs. Taking into account rather large random perturbations, the results show good robustness of the proposed method.

The 2D construction of circular arc and line approx from field road data are presented in Figs 6 and 7. In Fig. 6, one can find an example of an urban road with elements consisting of parts with small radii and short straight line segments, while Fig. 7 explains the meaning of additional information:

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

4. Conclusions

In this paper, the method for finding road centreline curvature from raw GPS data is presented. The approach consists of a few processing steps. First raw data of each road section using B-splines is fitted then the equidistant vertex of polyline of the fitted curve is generated. Using the least square method, the plane that best fits the points on the unit sphere, and find the circle that is the intersection of the plane and the unit sphere is found.

Stereographic projection of this circle back to the equatorial plane gives the corresponding circular arc and curvature.

The proposed method has been applied on a real 3D data model of the whole Slovenian state road network (more than 6000 km), and the performance is very satisfactory. Furthemore, the approach is simple and straightforward. Last but not least, it also works in higher dimensions, which is not the case with any other known methods. Estimation of clothoide transition lines still has to be explored.

doi: 10.3846/bjrbe.2011.21

Received 17 March 2010; accepted 2 September 2010

References

Barai, S. K. 2003. Data Mining Applications in Transportation Engineering, Transport 18(5): 216-223.

Bartels, R.; Beatty, J.; Barsky, B. 1987. An Introduction to Splines for Use in Computer Graphics and Geometric Modeling. Los Altos: Morgan Kaufmann, 476 p. ISBN 0934613273.

Cafiso, S.; Di Graziano, A.; Di Silvestro, G.; La Cava, G.; Persaud, B. 2010. Development of Comprehensive Accident Models for Two-Lane Rural Highways Using Exposure, Geometry, Consistency and Context Variables, Accident Analysis and Prevention 42(4): 1072-1079. doi:10.1016/j.aap.2009.12.015

Coeurjolly, D.; Gerard, J.; Reveilles, J. P.; Tougne, L. 2001. An Elementary Algorithm for Digital Arc Segmentation, Electronic Notes in Theoretical Computer Science 46: 355-370. doi:10.1016/S1571-0661(04)80997-6

Dell'Acqua, G.; Russo, F. 2010. Speed Factors on Low-Volume Roads for Horizontal Curves and Tangents, The Baltic Journal of .oad and Bridge Engineering 5(2): 89-97. doi:10.3846/bjrbe.2010.13

Discetti, P. 2010. Experimetal Analysis on Hairpin Curves, The Baltic Journal of Road and Bridge Engineering 5(3): 148-155. doi:10.3846/bjrbe.2010.21

Estrozi, L. F.; Campos, A. G.; Rios, L. G.; Cesar, R. M. 1999. Comparing Curvature Estimation Technique, in Proc. of the 4th Simposio Brasileiro de Automa?ao Inteligente (SBAI), Sao Paulo, Brazil, 58-63.

Fairney, D. P.; Fairney, P. T. 1994. On the Accuracy of Point Curvature Estimators in a Discrete Environment, Image and Vision Computing 12(5): 259-265. doi:10.1016/0262-8856(94)90031-0

Farin, G. 1997. Curves and Surfaces for Computer Aided Geometric Design: a Practical Guide. San Diego: Academic Press. 429 p. ISBN 0122490541.

Gontran, H.; Gillieron, P. Y.; Skaloud, J. 2005. Precise Road Geometry for Integrated Transport Safety Systems, in The 2nd Conference STRC 05, EPFL. Laboratoire de Topometrie.

Horng, J.-H.; Li, J. T. 2001. A Dynamic Programming Approach for Fitting Digital Planar Curves with Line Segments and Circular Arcs, Pattern Recognition letters 22(2): 183-197. doi:10.1016/S0167-8655(00)00104-5

Ichoku, C.; Deffontaines, B.; Chorowicz, J. 1996. Segmentation of Digital Plane Curves: a Dynamic Focusing Approach, Pattern Recognition Letters 17(7): 741-750. doi:10.1016/0167-8655(96)00015-3

Imran, M.; Hasan, Y.; Patterson, D. 2006. GPS-GIS Based Procedure for Tracking Vehicle Path on Horizontal Alignments, Computer-Aided Civil and Infrastructure Engineering 21(5): 383-394. doi:10.1111/j.1467-8667.2006.00444.x

Laurent, P. J.; Mehaute, A. L.; Schumaker, L. L. 1991. Curves and Surfaces. Boston: Academic Press. 514 p. ISBN 0124386601.

Lewiner, T.; Gomes, J. D.; Jr.; Lopes, H.; Craizer, M. 2004. ArcLength Based Curvature Estimator, in Proc. of the 17th Brazilian Symposium on Computer Graphics and Image. October 17-20, 2004. IEEE Press, 250-257. doi:10.1109/SIBGRA.2004.1352968

Mokhtarian, F.; Mackworth, A. 1992. A Theory of Multiscale, Curvature-Based Shape Representation for Planar Curves, IEEE Transactions on Pattern Analysis and Machine Intelligence 14(8): 789-805. doi:10.1109/34.149591

Pei, S. C.; Horng, J. H. 1996. Optimum Approximation of Digital Planar Curves Using Circular Arcs, Pattern Recognition 29(3): 383-388. doi:10.1016/0031-3203(95)00104-2

Pei, S. C.; Horng, J. H. 1995. Fitting Digital Curve Using Circular Arcs, Pattern Recognition 28(1): 107-116.

Perez, J. C.; Vidal, E. 1994. Optimum Polygonal Approximation of Digitized Curves, Pattern Recognition Letters 15(8): 743-750. doi:10.1016/0167-8655(94)90002-7

Pellegrino, O. 2009. An Analysis of the Effect of Roadway Design on Driver's Workload, The Baltic Journal of Road and Bridge Engineering 4(2): 45-53. doi:10.3846/1822-427X.2009.4.45-53

Perco, P. 2008. Influence of the General Character of Horizontal Alignment on Operating Speed of Two-Lane Rural Roads, Transportation Research Record 2075: 16-23. doi: 10.3141/2075-03

Pikaz, A.; Averbuch, A. 1996. On Automatic Threshold Selection for Polygonal Approximations of Digital Curves, Pattern Recognition 29(11): 1835-1845. doi:10.1016/0031-3203(96)00037-4

Ray, B. K.; Ray, K. S. 1994. A Non-Parametric Sequential Method for Polygonal Approximation of Digital Curves, Pattern Recognition Letters 15(2): 161-167. doi:10.1016/0167-8655(94)90045-0

Ray, B. K.; Ray, K. S. 1993. Determination of Optimal Polygon from Digital Curves Using L1 Norm, Pattern Recognition 26 (4): 505-509. doi:10.1016/0031-3203(93)90106-7

Rosin, P. L.; West, G. A. W. 1989. Segmentation of Edges into Lines and Arcs, Image and Vision Computing 7 (2): 109-114. doi:10.1016/0262-8856(89)90004-8

Schroedl, S.; Wagstaff, K.; Rogers, S.; Langley, P.; Wilson, C. 2004. Mining GPS Traces for Map Refinement, Data Mining and inowledge Discovery 9(1): 59-87. doi:10.1023/B:DAMI.0000026904.74892.89

Sliupas, T. 2009. The Impact of Road Parameters and the Surrounding Area on Traffic Accidents, Transport 24(1): 42-47. doi:10.3846/1648-4142.2009.24.42-47

Vorobjovas, V. 2011. Assurance of the Function of Low-Volume Roads for the Improvement of Driving Conditions, The Baltic Journal of Road and Bridge Engineering 6(1): 67-75. doi:10.3846/bjrbe.2011.09

Worring, M.; Smeulders, A. W. M. 1993. Digital Curvature Estimation, CVGIP: Image Understanding 58(3): 366-382. doi:10.1006/ciun.1993.1048

Yang, S. N.; Du, W. C. 1996. Numerical Methods for Approximating Digitized Curves by Piecewise Circular Arcs, Journal of Computational and Applied Mathematics 66(1-2): 557-569. doi:10.1016/0377-0427(95)00191-3

Zuriaga, A. M. P.; Garcia, A. G.; Torregrosa F. J. C.; D'Attoma, P. 2010. Modeling Operating Speed and Deceleration on Two-Lane Rural Roads with Global Positioning System Data, Transportation Research Record 2171: 11-20. doi: 10.3141/2171-02

Peter Lipar (1), Mitja Lakner (2), Tomaz Maher (3), Marijan Zura (4), ([mail])

Dept of Transport, University of Ljubljana, Jamova cesta 2, 1000 Ljubljana, Slovenia

E-mails: (1) plipar@fgg.uni-lj.si; (2) mlakner@fgg.uni-lj.si; (3) tmaher@fgg.uni-lj.si; (4) mzura@fgg.uni-lj.si

Printer friendly Cite/link Email Feedback | |

Author: | Lipar, Peter; Lakner, Mitja; Maher, Tomaz; Zura, Marijan |
---|---|

Publication: | The Baltic Journal of Road and Bridge Engineering |

Article Type: | Report |

Geographic Code: | 4EXSL |

Date: | Sep 1, 2011 |

Words: | 3490 |

Previous Article: | A proposed methodology for the management of low-volume roads in Spain/ Metodologija, siuloma mazo eismo intensyvumo keliams valdyti ispanijoje/... |

Next Article: | Development of Korean Pavement Design Guide for asphalt pavements based on the mechanistic-empirical design principle/Korejos asfalto dangu... |

Topics: |