Recognition and Reconstruction of Zebra Crossings on Roads from Mobile Laser Scanning Data

Lin Li 1,2,3, Da Zhang 1,4,*, Shen Ying 1,3,* and You Li 1 1 School of Resource and Environmental Sciences, Wuhan University, Wuhan 430079, China; lilin@whu.edu.cn (L.L.); boycecug@gmail.com (Y.L.) 2 Collaborative Innovation Center of Geospatial Technology, Wuhan University, Wuhan 430079, China 3 The Key Laboratory for Geographical Information System, Ministry of Education, Wuhan 430079, China 4 Power China Zhongnan Engineering Corporation Limited, Changsha 410014, China * Correspondence: zhangda_whu@163.com (D.Z.); shy@whu.edu.cn (S.Y.); Tel.: +86-131-6329-5128 (D.Z.)


Introduction
Road markings, as critical transportation infrastructure, provide drivers and pedestrians with information about traffic regulations, warnings, and guidance [1].The recognition and extraction of road markings have been seen as important functions in many fields, such as traffic safety management [2,3], driver assistance [4][5][6], and intelligent transportation [7,8].Traditional studies have mostly focused on digital images and videos [9][10][11][12][13][14][15].The extraction results are sometimes incomplete or insufficient due to poor weather conditions, lighting conditions, and complex shadowing from trees.In addition, the results fail to provide accurate three-dimensional (3D) coordinates of objects, which are crucial inputs to intelligent transportation systems and 3D city modelling.
Recent years have seen the emergence of mobile laser scanning (MLS) as a leading technology for extracting information about the surfaces of urban objects.MLS systems, which integrate laser scanners, global positioning system (GPS), inertial navigation system (INS), and charge-coupled device (CCD) cameras [16], collect information, such as 3D geospatial, texture, and laser intensity data, from complex urban areas when a vehicle is on the move.Such systems have become a promising and cost-effective solution for rapid road environment modelling.Most methods proposed in previous studies are designed for point cloud classification [17][18][19][20], building footprint extraction, façade reconstruction [21][22][23], and detection of vertical pole-like objects [24][25][26] in a road environment.Only a few studies have explored the recognition and extraction of road markings.Jaakkola et al. [27] first generated georeferenced feature images according to elevation and intensity by applying an interpolation method and then segmenting road markings and curbstones by applying thresholding and morphological operations to elevation and intensity images.Yang et al. [28] calculated the weights and pixel values of grey images based on spatial distributions (e.g., planar distance, elevation difference, and point density) of laser scanning points, which improved their algorithm for generating feature images.Then, they applied an intensity filter and an elevation filter, followed by constraints on shape and distribution patterns.The above-mentioned methods transform 3D points into 2D images because addressing a large volume of unorganized points is time-consuming and complex.The transformation improves the computational efficiency and enables one to capitalize on well-established image processing methods.However, this transformation also causes roughness in detail, especially when extracting small objects, such as road markings.
Kammel [29] and Chen [30] applied the Radon transform and Hough transform to extract solid-edge-line and dashed-lane-line markings, respectively, from MLS points.These methods are effective when extracting straight markings; however, they exhibit a weakness in extracting curve markings.Since curve markings are usually irregular, it is difficult to choose a suitable curve model.
Simple global intensity thresholding is often used to extract road markings [27,28,31]; however, the markings' non-uniform intensity makes this method less effective in some cases because intensity values are affected by material, laser incidence angle, and range.Guan et al. [32] proposed a novel method that segments the intensity images with multiple thresholds related to the point density.Using their method, an image was partitioned into several blocks in accordance with the point density distribution characteristics.Within different blocks, local optimal thresholds were estimated to extract road markings.However, substantial noise was introduced in this method.
In addition to the extraction of road markings, the recognition of types is also a necessary and challenging task, especially for zebra crossings, which are located in urban road intersections and have important functions in traffic safety management.Mancini et al. [33] identified zebra crossings with area, perimeter, and length-width ratios following connected component labelling.Riveiro et al. [34] used a Canny edge detector and standard Hough transform to detect a set of parallel lines that have similar directions as the road centreline.Yu et al. [35] distinguished zebra crossings from other rectangular-shaped markings according to the geometric perpendicularity of their distribution directions and road centrelines.These studies mostly focused on detecting stripes and did not provide specific information about the areas.For guide systems for blind individuals, mobile robot navigation, etc., it is impossible to confirm whether the frontal area is a zebra crossing without such information.Another problem is that the method is invalid when the distribution directions of zebra crossings and road centrelines are not vertical.
To overcome the aforementioned limitations, we propose a stepwise procedure for recognizing and reconstructing zebra crossings using mobile laser scanning data in this paper.The contributions of this paper are as follows: (1) an adaptive thresholding method based on road surface partitioning was designed to compensate for non-uniformities in intensity data and extract all types of road markings; (2) a dispersion degree filtering method was applied to reduce the noise; and (3) zebra crossings are recognized and reconstructed according to geometrical features, so that we obtain more specific information about the area, including start positions, end positions, distribution directions of zebra crossings, and road centreline directions.
The remainder of this paper is organized as follows: the stepwise description of the proposed method is presented in Section 2; in Section 3, we test the proposed method on MLS data captured in Wuhan, China; following the experiments, we discuss the results; and, finally, conclusions are drawn in Section 4.

Method
The method consists of three main procedures: road surface segmentation, road marking extraction, and zebra crossing recognition and reconstruction.Figure 1  The method consists of three main procedures: road surface segmentation, road marking extraction, and zebra crossing recognition and reconstruction.Figure 1 illustrates the complete experimental procedures used in this study.

Road Surface Segmentation
In an urban environment, road surfaces are generally flat, with small elevation jumps caused by curbstones on the road boundaries, as shown in Figure 2. The elevations of road boundary points change substantially more quickly than do those of road surface points.The gradient of a scalar field reflects the rate of change of the scalar; therefore, we attempt to separate roads from other points via elevation gradients and apply a region growing algorithm to the elevation-gradient feature image for precise road surface segmentation.

Road Surface Segmentation
In an urban environment, road surfaces are generally flat, with small elevation jumps caused by curbstones on the road boundaries, as shown in Figure 2. The elevations of road boundary points change substantially more quickly than do those of road surface points.The gradient of a scalar field reflects the rate of change of the scalar; therefore, we attempt to separate roads from other points via elevation gradients and apply a region growing algorithm to the elevation-gradient feature image for precise road surface segmentation.
ISPRS Int.J. Geo-Inf.2016, 5, 125 3 of 16 The method consists of three main procedures: road surface segmentation, road marking extraction, and zebra crossing recognition and reconstruction.Figure 1 illustrates the complete experimental procedures used in this study.

Road Surface Segmentation
In an urban environment, road surfaces are generally flat, with small elevation jumps caused by curbstones on the road boundaries, as shown in Figure 2. The elevations of road boundary points change substantially more quickly than do those of road surface points.The gradient of a scalar field reflects the rate of change of the scalar; therefore, we attempt to separate roads from other points via elevation gradients and apply a region growing algorithm to the elevation-gradient feature image for precise road surface segmentation.

Elevation Filtering
Trees next to roads, street lamps, etc., may cause difficulties in extracting accurate road surface points.For instance, tree crowns can cover road boundaries when laser scanning points are projected onto the XY plane.Therefore, we extend the histogram concavity analysis algorithm proposed by Rosenfeld [36] and use the extended algorithm to select an appropriate threshold for elevation filtering.The algorithm is applicable to both unimodal and bimodal histograms.
As shown in Figure 4, we first draw a histogram hs based on F, one of the laser scanning points' feature properties.The class width is w, and each rectangle is numbered i = 1 ,..., n.We added two points, (f1, 0) and (fn, 0), for defining region hs more conveniently and accurately.For a rectangle i, let the feature value fi be the X coordinate of the upper side's midpoint, and let hi be the height: To find the concavities of hs, we first construct its convex hull HS.This is the smallest convex polygon containing (f1, 0), (fn, 0) and (fi, hi) (i = 1,2, …, n), midpoints of all the rectangle's upper sides.Hi is the height of HS when the feature value is fi.The depth of the concavity Di is determined as follows: Any concavity in the histogram may be the location of a threshold; however, not all points of the concavity are good candidates.The deeper the concavity depth is, the larger the differences in objects' feature values between both sides; thus, we shall consider those points for which Di is a local maximum as candidates for positions of the threshold T:

Elevation Filtering
Trees next to roads, street lamps, etc., may cause difficulties in extracting accurate road surface points.For instance, tree crowns can cover road boundaries when laser scanning points are projected onto the XY plane.Therefore, we extend the histogram concavity analysis algorithm proposed by Rosenfeld [36] and use the extended algorithm to select an appropriate threshold for elevation filtering.The algorithm is applicable to both unimodal and bimodal histograms.
As shown in Figure 4, we first draw a histogram hs based on F, one of the laser scanning points' feature properties.The class width is w, and each rectangle is numbered i = 1 ,..., n.We added two points, (f 1 , 0) and (f n , 0), for defining region hs more conveniently and accurately.

Elevation Filtering
Trees next to roads, street lamps, etc., may cause difficulties in extracting accurate road surface points.For instance, tree crowns can cover road boundaries when laser scanning points are projected onto the XY plane.Therefore, we extend the histogram concavity analysis algorithm proposed by Rosenfeld [36] and use the extended algorithm to select an appropriate threshold for elevation filtering.The algorithm is applicable to both unimodal and bimodal histograms.
As shown in Figure 4, we first draw a histogram hs based on F, one of the laser scanning points' feature properties.The class width is w, and each rectangle is numbered i = 1 ,..., n.We added two points, (f1, 0) and (fn, 0), for defining region hs more conveniently and accurately.For a rectangle i, let the feature value fi be the X coordinate of the upper side's midpoint, and let hi be the height: To find the concavities of hs, we first construct its convex hull HS.This is the smallest convex polygon containing (f1, 0), (fn, 0) and (fi, hi) (i = 1,2, …, n), midpoints of all the rectangle's upper sides.Hi is the height of HS when the feature value is fi.The depth of the concavity Di is determined as follows: Any concavity in the histogram may be the location of a threshold; however, not all points of the concavity are good candidates.The deeper the concavity depth is, the larger the differences in objects' feature values between both sides; thus, we shall consider those points for which Di is a local maximum as candidates for positions of the threshold T: For a rectangle i, let the feature value f i be the X coordinate of the upper side's midpoint, and let h i be the height: To find the concavities of hs, we first construct its convex hull HS.This is the smallest convex polygon containing (f 1 , 0), (f n , 0) and (f i , h i ) (i = 1, 2, . . ., n), midpoints of all the rectangle's upper sides.H i is the height of HS when the feature value is f i .The depth of the concavity D i is determined as follows: Any concavity in the histogram may be the location of a threshold; however, not all points of the concavity are good candidates.The deeper the concavity depth is, the larger the differences in objects' feature values between both sides; thus, we shall consider those points for which D i is a local maximum as candidates for positions of the threshold T: We select elevation information of point clouds as feature properties for histogram concavity analysis in this study.In urban road environments, there are usually a large number of road surface points, and their distribution is concentrated.Trees, street lamps, and other objects that cause serious interference in road segmentation have greater heights, and the points are dispersed.Therefore, based on the elevation distribution characteristics of point clouds, we conclude that the peak of the histogram corresponds to the road surface and that "shoulders" on the right side of the peak correspond to those objects that could cause interference in road segmentation.As shown in Figure 5c, the possible threshold locates in the "shoulder" on the right.First, let h f be the elevation corresponding to the largest frequency, and let h max be the highest elevation of all points.Then, calculate the threshold T in the range [h f , h max ] by analysing the histogram concavity.Finally, we filter out points with an elevation of less than T and use the remaining points for subsequent procedures.
ISPRS Int.J. Geo-Inf.2016, 5, 125 5 of 16 We select elevation information of point clouds as feature properties for histogram concavity analysis in this study.In urban road environments, there are usually a large number of road surface points, and their distribution is concentrated.Trees, street lamps, and other objects that cause serious interference in road segmentation have greater heights, and the points are dispersed.Therefore, based on the elevation distribution characteristics of point clouds, we conclude that the peak of the histogram corresponds to the road surface and that "shoulders" on the right side of the peak correspond to those objects that could cause interference in road segmentation.As shown in Figure 5(c), the possible threshold locates in the "shoulder" on the right.First, let hf be the elevation corresponding to the largest frequency, and let hmax be the highest elevation of all points.Then, calculate the threshold T in the range [hf, hmax] by analysing the histogram concavity.Finally, we filter out points with an elevation of less than T and use the remaining points for subsequent procedures.

Segmentation by Region-Growing
To improve the computation speed of the proposed method and apply image processing algorithms, laser scanning points are projected onto an XY plane to generate a georeferenced feature image I.The grey value of each cell equals the central point's elevation calculated using the Inverse Distance Weighted (IDW) interpolation method [37] and normalization.Then, the bilinear interpolation method is used to smooth the image because there are no points in some cells, which causes image noise.Finally, the elevation feature image I is converted into an elevation-gradient feature image G as follows:

Segmentation by Region-Growing
To improve the computation speed of the proposed method and apply image processing algorithms, laser scanning points are projected onto an XY plane to generate a georeferenced feature image I.The grey value of each cell equals the central point's elevation calculated using the Inverse Distance Weighted (IDW) interpolation method [37] and normalization.Then, the bilinear interpolation method is used to smooth the image because there are no points in some cells, which causes image noise.Finally, the elevation feature image I is converted into an elevation-gradient feature image G as follows: Since the elevation gradients of the road surface and road boundary are sufficiently distinct, we use a single threshold T G for image binarization.In the binarized image B, the pixel value is set to 1 if the corresponding pixel value in G is greater than T G ; otherwise, it is set to 0. The result is shown as Figure 6b.For bridging gaps on road boundaries, dilation is applied to binarized images.Figure 6c shows the results of dilating the image B with a 3 ˆ3 structuring element, and the road surface is the black area surrounded by a clear white road boundary.if the corresponding pixel value in G is greater than TG; otherwise, it is set to 0. The result is shown as Figure 6b.For bridging gaps on road boundaries, dilation is applied to binarized images.Figure 6c shows the results of dilating the image B with a 3 × 3 structuring element, and the road surface is the black area surrounded by a clear white road boundary.Region-growing is an effective method for extracting the road surface in an image.The first step is to choose a point in the trajectory data and confirm that the value of a pixel in which the point locates is 0.Then, the pixel can be determined as a seed point.All of the pixels that are 8-connected to the seed point and whose values are also 0 are appended to the seed point to form a larger region.8-connected pixels are neighbours to every pixel that touches one of their edges or corners.Then, we dilate the result of the region-growing process using a 3 × 3 structuring element to compensate for the error caused by the previous dilation.Finally, the point clouds of road surfaces are converted based on the optimized region-growing results.

Adaptive Thresholding Based on Road Surface Partitioning
Usually, a road surface is composed of asphalt and concrete and exhibits low or diffuse reflection properties when subject to an incident laser.Road markings are highly reflective white or yellow coatings painted on the road surface.Objects with higher reflectances correspond to stronger laser signals; therefore, the laser intensity value is a key feature for distinguishing road markings from road surfaces.However, the intensity values are also affected by the laser incidence angles and the distances between the target and the scanner centre, which makes single global thresholding less effective for segmentation.Therefore, adaptive intensity thresholding based on road surface partitioning is proposed to solve the problems caused by non-uniform intensities.
Generally, the farther away from the trajectory data, the lower the intensity of road markings.The materials that constitute different road sections also vary considerably.Thus, the road surface is partitioned into non-overlapping rectangles Recti, as shown in Figure 7.The X axis is the direction of the vehicle trajectory, and the Y axis is perpendicular to the X axis in the horizontal plane.The length and width of the rectangle are Rx and Ry, respectively.The size of the rectangle is related to the extent of intensity evenness.Rx and Ry should be set smaller when the intensity distributes more unevenly in order to ensure a uniform intensity distribution in each rectangle.Region-growing is an effective method for extracting the road surface in an image.The first step is to choose a point in the trajectory data and confirm that the value of a pixel in which the point locates is 0.Then, the pixel can be determined as a seed point.All of the pixels that are 8-connected to the seed point and whose values are also 0 are appended to the seed point to form a larger region.8-connected pixels are neighbours to every pixel that touches one of their edges or corners.Then, we dilate the result of the region-growing process using a 3 ˆ3 structuring element to compensate for the error caused by the previous dilation.Finally, the point clouds of road surfaces are converted based on the optimized region-growing results.

Adaptive Thresholding Based on Road Surface Partitioning
Usually, a road surface is composed of asphalt and concrete and exhibits low or diffuse reflection properties when subject to an incident laser.Road markings are highly reflective white or yellow coatings painted on the road surface.Objects with higher reflectances correspond to stronger laser signals; therefore, the laser intensity value is a key feature for distinguishing road markings from road surfaces.However, the intensity values are also affected by the laser incidence angles and the distances between the target and the scanner centre, which makes single global thresholding less effective for segmentation.Therefore, adaptive intensity thresholding based on road surface partitioning is proposed to solve the problems caused by non-uniform intensities.
Generally, the farther away from the trajectory data, the lower the intensity of road markings.The materials that constitute different road sections also vary considerably.Thus, the road surface is partitioned into non-overlapping rectangles Rect i , as shown in Figure 7.The X axis is the direction of the vehicle trajectory, and the Y axis is perpendicular to the X axis in the horizontal plane.The length and width of the rectangle are R x and R y , respectively.The size of the rectangle is related to the extent of intensity evenness.R x and R y should be set smaller when the intensity distributes more unevenly in order to ensure a uniform intensity distribution in each rectangle.There are two possibilities concerning the number of point types in a rectangle: (a) only one type, i.e., road surface points, or (b) two types, i.e., road surface points and road marking points.Otsu's algorithm [38] is first used to find the optimal intensity threshold in each rectangle.Then, the two cases are separated based on the thresholding results.
The point set PA = { p1, p2, …, pm } represents the points whose intensities are larger than the threshold, and PB = { p1, p2, …, pm } represents the remaining points in the rectangle.The 3D coordinates of pi are (xi, yi, zi), and its intensity value is Ii.
For case (a), the points in PA and PB are both road surface points.For case (b), the points in PA are road marking points, and the points in PB are road surface points.The distance of cluster centres between PA and PB in case (b) is far larger than that in case (a).This distance dI is calculated as follows: The ratio of the number of points between PA and PB is also a critical element to the judgement of cases.In Figure 7, Rect1 and Rect2 are examples of case (a) and case (b).Their thresholding results are shown in Figure 8.In case (b), the number of points in PA is substantially larger than that in PB, which leads to a high ratio.The ratio is defined as follows: According to the above analysis, the two cases can be distinguished using the following formula: where Td and Tr are the thresholds of the cluster-centre distance and the ratio of the number of points, respectively.There are two possibilities concerning the number of point types in a rectangle: (a) only one type, i.e., road surface points, or (b) two types, i.e., road surface points and road marking points.Otsu's algorithm [38] is first used to find the optimal intensity threshold in each rectangle.Then, the two cases are separated based on the thresholding results.
The point set P A = { p 1 , p 2 , . . ., p m } represents the points whose intensities are larger than the threshold, and P B = { p 1 , p 2 , . . ., p m } represents the remaining points in the rectangle.The 3D coordinates of p i are (x i , y i , z i ), and its intensity value is I i .
For case (a), the points in P A and P B are both road surface points.For case (b), the points in P A are road marking points, and the points in P B are road surface points.The distance of cluster centres between P A and P B in case (b) is far larger than that in case (a).This distance d I is calculated as follows: The ratio of the number of points between P A and P B is also a critical element to the judgement of cases.In Figure 7, Rect 1 and Rect 2 are examples of case (a) and case (b).Their thresholding results are shown in Figure 8.In case (b), the number of points in P A is substantially larger than that in P B , which leads to a high ratio.The ratio is defined as follows: According to the above analysis, the two cases can be distinguished using the following formula: @Rect i : # i f pd I ă T d q and pratio ą T r q , case paq otherwise, case pbq where T d and T r are the thresholds of the cluster-centre distance and the ratio of the number of points, respectively.Finally, for case (b), all of the points in PA are reserved as the coarse results obtained from road marking extraction.

Dispersion Degree Filtering
Some parts of road surfaces have similar material properties as those of road markings, which causes noise in the coarse extraction results, as shown in Figure 9.The road marking points are more concentrated than the noise; therefore, noise is proposed to be removed according to the difference in dispersion degrees.The dispersion degree Dp of a point p(x, y, z) is defined as follows: (8) where Np denotes the number of local neighbourhood points.By removing the points whose dispersion degrees are larger than the threshold TD, accurately extracted road markings can be obtained.

The Model of Zebra Crossings
A zebra crossing is an area that consists of a group of broad white stripes painted on the road.As shown in Figure 10, our model of a zebra crossing contains the following four elements: L1 and L2 define the start and end positions, respectively; Vr is the direction of the road centreline, which is also the direction along which vehicles travel; and Vz is the distribution direction of the zebra crossing, which guides pedestrians to cross the road safely.Finally, for case (b), all of the points in P A are reserved as the coarse results obtained from road marking extraction.

Dispersion Degree Filtering
Some parts of road surfaces have similar material properties as those of road markings, which causes noise in the coarse extraction results, as shown in Figure 9.The road marking points are more concentrated than the noise; therefore, noise is proposed to be removed according to the difference in dispersion degrees.Finally, for case (b), all of the points in PA are reserved as the coarse results obtained from road marking extraction.

Dispersion Degree Filtering
Some parts of road surfaces have similar material properties as those of road markings, which causes noise in the coarse extraction results, as shown in Figure 9.The road marking points are more concentrated than the noise; therefore, noise is proposed to be removed according to the difference in dispersion degrees.The dispersion degree Dp of a point p(x, y, z) is defined as follows: (8) where Np denotes the number of local neighbourhood points.By removing the points whose dispersion degrees are larger than the threshold TD, accurately extracted road markings can be obtained.

The Model of Zebra Crossings
A zebra crossing is an area that consists of a group of broad white stripes painted on the road.As shown in Figure 10, our model of a zebra crossing contains the following four elements: L1 and L2 define the start and end positions, respectively; Vr is the direction of the road centreline, which is also the direction along which vehicles travel; and Vz is the distribution direction of the zebra crossing, which guides pedestrians to cross the road safely.The dispersion degree D p of a point p(x, y, z) is defined as follows: where N p denotes the number of local neighbourhood points.By removing the points whose dispersion degrees are larger than the threshold T D , accurately extracted road markings can be obtained.

The Model of Zebra Crossings
A zebra crossing is an area that consists of a group of broad white stripes painted on the road.As shown in Figure 10, our model of a zebra crossing contains the following four elements: L 1 and L 2 define the start and end positions, respectively; V r is the direction of the road centreline, which is also the direction along which vehicles travel; and V z is the distribution direction of the zebra crossing, which guides pedestrians to cross the road safely.

Detection of Zebra Stripes
Design standards in most countries provide regulations on the exact size and shape of road markings.It is advantageous to recognize different types of road markings based on their different sizes and shapes.The stripe is actually a rectangle with a fixed size.Therefore, there are two factors we can use to recognize stripes: (a) rectangular features and (b) fixed lengths Lz and widths Wz.
To cluster neighbouring road marking points, intensity feature images transformed from point clouds are binarized, followed by 8-connected component labelling.For each connected region, an ellipse with the same second moments as the 8-connected region is calculated.e is the eccentricity of the ellipse, and its value ranges from 0 to 1.The closer this value is to 1, the more likely the region is a rectangle.Based on experience, regions whose e values are larger than 0.99 are candidate stripes, and the corresponding points are denoted as the point set P = {p1, p2, …, pk}.
To calculate the length and width of the candidate stripes, the principal component analysis (PCA) method is performed to judge the principal distribution direction of points in P on the XY plane.
For all the points in P, the correlation between xi and yi can be determined through their variances as follows: where x and y are the average values of xi and yi.
The covariance matrix C can be established using Equation (10): As shown in Figure 11, through the eigenvalue decomposition (EVD) of C, the eigenvector associated with the larger eigenvalue V1 is defined as the first principal direction, which is parallel to the long side of rectangle; the eigenvector associated with the smaller eigenvalue V2 is defined as the second principal direction, which is parallel to the short side of the rectangle.One dimension matrixes M1 and M2 are established by projecting all points' plane coordinates to V1 and V2, respectively.Then, the length Lz and width Wz of the region can be calculated as follows: Figure 10.The model of zebra crossings: (a) the distribution direction is perpendicular to the road centreline; and (b) the distribution direction is oblique to the road centreline.

Detection of Zebra Stripes
Design standards in most countries provide regulations on the exact size and shape of road markings.It is advantageous to recognize different types of road markings based on their different sizes and shapes.The stripe is actually a rectangle with a fixed size.Therefore, there are two factors we can use to recognize stripes: (a) rectangular features and (b) fixed lengths L z and widths W z .
To cluster neighbouring road marking points, intensity feature images transformed from point clouds are binarized, followed by 8-connected component labelling.For each connected region, an ellipse with the same second moments as the 8-connected region is calculated.e is the eccentricity of the ellipse, and its value ranges from 0 to 1.The closer this value is to 1, the more likely the region is a rectangle.Based on experience, regions whose e values are larger than 0.99 are candidate stripes, and the corresponding points are denoted as the point set P = {p 1 , p 2 , . . ., p k }.
To calculate the length and width of the candidate stripes, the principal component analysis (PCA) method is performed to judge the principal distribution direction of points in P on the XY plane.
For all the points in P, the correlation between x i and y i can be determined through their variances as follows: cov px, yq " where x and y are the average values of x i and y i .
The covariance matrix C can be established using Equation (10): As shown in Figure 11, through the eigenvalue decomposition (EVD) of C, the eigenvector associated with the larger eigenvalue V 1 is defined as the first principal direction, which is parallel to the long side of rectangle; the eigenvector associated with the smaller eigenvalue V 2 is defined as the second principal direction, which is parallel to the short side of the rectangle.One dimension matrixes M 1 and M 2 are established by projecting all points' plane coordinates to V 1 and V 2 , respectively.Then, the length L z and width W z of the region can be calculated as follows:

#
L z " max pM 1 q ´min pM 1 q W z " max pM 2 q ´min pM 2 q (11) The regions whose length and width satisfy the design standards can be reserved as zebra stripes.Considering wear on road markings and calculation errors in real-world cases, the value ranges of Lz and Wz are set as follows: where Ls and Ws are the standard length and width of zebra stripes.

Reconstruction of Zebra Crossings
The centroids of all stripes in a zebra crossing locate along a straight line.This line is found on the centre axis of the zebra crossing, which is important to area reconstruction.
Random sample consensus (RANSAC) [39] is an effective iterative algorithm for mathematical model fitting, such as linear fitting.By adjusting the number of iterations nR and the residual threshold TR, optimal model parameters can be estimated from a set of observed data containing noise.This is adopted to solve the problem of the fitting of zebra crossing's centre axes in this paper.
We first calculate the coordinates of all stripe centroids on the XY plane.Then, the RANSAC algorithm is applied to these centroids.To ensure the accuracy of the results, there should be at least three points used for fitting with the estimated linear model.Finally, we directly obtain some important information: (a) the number of zebra crossings: the number of iterations; (b) the centre axis of zebra crossings: the lines fitted by RANSAC; and (c) stripes belonging to the same zebra crossing.The distribution direction Vz is the same as the direction of the centre axis.The direction of the road's centreline Vr is calculated by averaging the stripes' first principle direction in a zebra crossing.L1 and L2 are obtained by translating the centre axis along Vr, where the translational distance is Ls/2.This completes the recognition and reconstruction of zebra crossings.

Results and Discussion
The point clouds used in the experiment were captured by an Optech Lynx mobile mapping system, which consists of two laser scanners, one GPS receiver, and an inertial measurement unit.The original data is given in the WGS-84 coordinate system.Then the data is transformed from longitude and latitude coordinates to mathematical X and Y planar coordinates using Gauss projection.The survey area is in Guanggu, a part of the City of Wuhan, which is a major city in central China.
Figure 12 shows the three datasets selected for the evaluation of the performance of the proposed method.The figure contains vegetation (e.g., trees and bushes), street lamps, power lines, and cars in these areas.The roads in these datasets consist of straight sections, curved sections, and crossroads.Detailed information, including road length and the number of points, is presented in Table 1.The regions whose length and width satisfy the design standards can be reserved as zebra stripes.Considering wear on road markings and calculation errors in real-world cases, the value ranges of L z and W z are set as follows: where L s and W s are the standard length and width of zebra stripes.

Reconstruction of Zebra Crossings
The centroids of all stripes in a zebra crossing locate along a straight line.This line is found on the centre axis of the zebra crossing, which is important to area reconstruction.
Random sample consensus (RANSAC) [39] is an effective iterative algorithm for mathematical model fitting, such as linear fitting.By adjusting the number of iterations n R and the residual threshold T R , optimal model parameters can be estimated from a set of observed data containing noise.This is adopted to solve the problem of the fitting of zebra crossing's centre axes in this paper.
We first calculate the coordinates of all stripe centroids on the XY plane.Then, the RANSAC algorithm is applied to these centroids.To ensure the accuracy of the results, there should be at least three points used for fitting with the estimated linear model.Finally, we directly obtain some important information: (a) the number of zebra crossings: the number of iterations; (b) the centre axis of zebra crossings: the lines fitted by RANSAC; and (c) stripes belonging to the same zebra crossing.The distribution direction V z is the same as the direction of the centre axis.The direction of the road's centreline V r is calculated by averaging the stripes' first principle direction in a zebra crossing.L 1 and L 2 are obtained by translating the centre axis along V r , where the translational distance is ˘Ls /2.This completes the recognition and reconstruction of zebra crossings.

Results and Discussion
The point clouds used in the experiment were captured by an Optech Lynx mobile mapping system, which consists of two laser scanners, one GPS receiver, and an inertial measurement unit.The original data is given in the WGS-84 coordinate system.Then the data is transformed from longitude and latitude coordinates to mathematical X and Y planar coordinates using Gauss projection.The survey area is in Guanggu, a part of the City of Wuhan, which is a major city in central China.
Figure 12 shows the three datasets selected for the evaluation of the performance of the proposed method.The figure contains vegetation (e.g., trees and bushes), street lamps, power lines, and cars in these areas.The roads in these datasets consist of straight sections, curved sections, and crossroads.Detailed information, including road length and the number of points, is presented in Table 1.

Segmentation of Road Surfaces
To section the experimental data into a number of blocks, we chose d=50m in dataset 2. In the other two datasets, we used d=30m because there are more curves and ramps.For each block histogram concavity analysis was used to obtain elevation thresholds.Then, following elevation filtering, point clouds were converted into elevation-gradient feature images.The grid size is a critical parameter in image generation.When the size is too small, only a few points, or, possibly, no points, fall inside the grids, whereas a large size may result in low image quality.Taking dataset 1 as an example, a block of data was selected to generate elevation-gradient images with different grid sizes of 0.05, 0.07, 0.09, and 0.11 m. Figure 13 presents the comparison results.Visual inspection suggests that there are few noise points on the road surface and that the details are clear when the grid size is 0.09 m; therefore, this value was applied in the experiment.The grid sizes used in dataset 2 and 3 were set to 0.12 m and 0.10 m, respectively, in the same way.

Segmentation of Road Surfaces
To section the experimental data into a number of blocks, we chose d = 50 m in dataset 2. In the other two datasets, we used d = 30 m because there are more curves and ramps.For each block histogram concavity analysis was used to obtain elevation thresholds.Then, following elevation filtering, point clouds were converted into elevation-gradient feature images.The grid size is a critical parameter in image generation.When the size is too small, only a few points, or, possibly, no points, fall inside the grids, whereas a large size may result in low image quality.Taking dataset 1 as an example, a block of data was selected to generate elevation-gradient images with different grid sizes of 0.05, 0.07, 0.09, and 0.11 m. Figure 13 presents the comparison results.Visual inspection suggests that there are few noise points on the road surface and that the details are clear when the grid size is 0.09 m; therefore, this value was applied in the experiment.The grid sizes used in dataset 2 and 3 were set to 0.12 m and 0.10 m, respectively, in the same way.

Segmentation of Road Surfaces
To section the experimental data into a number of blocks, we chose d=50m in dataset 2. In the other two datasets, we used d=30m because there are more curves and ramps.For each block histogram concavity analysis was used to obtain elevation thresholds.Then, following elevation filtering, point clouds were converted into elevation-gradient feature images.The grid size is a critical parameter in image generation.When the size is too small, only a few points, or, possibly, no points, fall inside the grids, whereas a large size may result in low image quality.Taking dataset 1 as an example, a block of data was selected to generate elevation-gradient images with different grid sizes of 0.05, 0.07, 0.09, and 0.11 m. Figure 13 presents the comparison results.Visual inspection suggests that there are few noise points on the road surface and that the details are clear when the grid size is 0.09 m; therefore, this value was applied in the experiment.The grid sizes used in dataset 2 and 3 were set to 0.12 m and 0.10 m, respectively, in the same way.For the binarization of elevation-gradient images, a threshold should be determined.The grey values of road surfaces generally range from 0 to 0.005, and the grey values of road boundaries are approximately 0.015; therefore, any value between 0.005 and 0.015 could be set as the threshold and we selected 0.005.
Finally, we segmented the road surfaces using a region-growing method; then, the 3D points associated with road surfaces could be extracted easily, as shown in Figure 14.The close-up views in the black rectangles indicate that the road surfaces are basically extracted accurately and completely.

Extraction of Road Markings
Several parameters and their values used in the extraction of road markings are listed in Table 2.They were mainly selected through a set of tests or based on prior knowledge.Then, road markings were extracted directly with adaptive thresholding and dispersion degree filtering from road surface points.All types of road markings could be extracted fairly well.However, a few road markings were abraded by cars and pedestrians, leading to some of the extraction results being incomplete.Figure 15 shows a part of the extracted road markings, including solid lines, dotted lines, arrow markings, and diamond markings.

Extraction of Road Markings
Several parameters and their values used in the extraction of road markings are listed in Table 2.They were mainly selected through a set of tests or based on prior knowledge.Then, road markings were extracted directly with adaptive thresholding and dispersion degree filtering from road surface points.All types of road markings could be extracted fairly well.However, a few road markings were abraded by cars and pedestrians, leading to some of the extraction results being incomplete.Figure 15 shows a part of the extracted road markings, including solid lines, dotted lines, arrow markings, and diamond markings.

Recognition and Reconstruction of Zebra Crossings
The standard length Ls and width Ws of the zebra stripes are 6 m and 0.4 m, respectively, in the three datasets, which satisfy the design standards of zebra crossings in China.After recognizing stripes according to the above standards, the RANSAC algorithm was applied to the centroids of

Recognition and Reconstruction of Zebra Crossings
The standard length L s and width W s of the zebra stripes are 6 m and 0.4 m, respectively, in the three datasets, which satisfy the design standards of zebra crossings in China.After recognizing stripes according to the above standards, the RANSAC algorithm was applied to the centroids of stripes with an n R of 5000 and a T R of 0.25 to obtain comprehensive information about zebra crossing areas.
A comparative study was conducted to compare our proposed zebra crossing recognition method with a recently published method: Riveiro's method [34].As listed in Table 3, a total of eleven zebra crossings were recognized at a recognition rate of 90.91% with method, which outperforms the other method.One zebra crossing was not detected due to the low reflectivity of road markings caused by serious abrasion, which decreases the completeness of road marking extraction, as shown in the Figure 16.

Recognition and Reconstruction of Zebra Crossings
The standard length Ls and width Ws of the zebra stripes are 6 m and 0.4 m, respectively, in the three datasets, which satisfy the design standards of zebra crossings in China.After recognizing stripes according to the above standards, the RANSAC algorithm was applied to the centroids of stripes with an nR of 5000 and a TR of 0.25 to obtain comprehensive information about zebra crossing areas.
A comparative study was conducted to compare our proposed zebra crossing recognition method with a recently published method: Riveiro's method [34].As listed in Table 3, a total of eleven zebra crossings were recognized at a recognition rate of 90.91% with our method, which outperforms the other method.One zebra crossing was not detected due to the low reflectivity of road markings caused by serious abrasion, which decreases the completeness of road marking extraction, as shown in the Figure 16.To further quantitatively evaluate the performance of our method, four measures were computed for each zebra crossing based on manually-extracted results.z and r represent the angle To further quantitatively evaluate the performance of our method, four measures were computed for each zebra crossing based on manually-extracted results.θ z and θ r represent the angle deviations of a zebra crossing's distribution direction and a road centreline's direction, respectively.The completeness r is used to describe how complete the detected zebra crossing areas are, and the correctness p is used to indicate what percentage of the detected zebra crossing areas are valid.r and p are defined as follows: r " TP{AP p " TP{VP where TP, AP, and VP are the number of road surface points belonging to (1) the correctly detected zebra crossing areas using the proposed method; (2) the zebra crossing areas collected using manual visual interpretation; and (3) the whole detected zebra crossing areas using the proposed method, respectively.
As shown by the quality evaluation results in Table 4, the completeness and correctness of recognized zebra crossings are both greater than 90%, the value of θ z is no greater than 2.5 ˝, and the maximum value of θ r is 1.2 ˝.In summary, our proposed method exhibits good performance in recognizing and reconstructing zebra crossings.

Conclusions
In this paper, we have proposed an effective method for recognizing and reconstructing zebra crossings using mobile laser scanning data.The proposed method first converts point clouds into elevation-gradient images and subsequently applies region-growing-based road surface segmentation.Second, road marking points are extracted with adaptive intensity thresholding based on road surface partitioning and dispersion degree filtering.Finally, the zebra crossing areas are recognized and reconstructed according to geometrical features.
The three datasets acquired by an Optech Lynx mobile mapping system were used to validate our zebra crossing recognition and reconstruction method.The experimental results demonstrate that the proposed method performs well and obtains high completeness and correctness values.The experiment has indicated three main advantages of the method: (1) the method is effective even when points of zebra crossings are incomplete; (2) the method can be effective when the distribution direction of zebra crossings and road centrelines are at arbitrary angles; and (3) more comprehensive information about zebra crossing areas, such as the extent of the area, is obtained.
These research findings could contribute to a more rapid, cost-effective, and comprehensive approach to traffic management and ensure maximum safety conditions for road users.However, our method could only be applied for post-processing instead of real-time use at present, because some parameters need to be selected based on prior knowledge or a set of tests.In the future, we will make further study on the algorithm of selecting optimal parameters automatically.Additionally, it is also important to enhance the computing efficiency of our method, because point clouds with better resolution and higher density are needed if we want to obtain more detailed information about urban objects.

Figure 1 .
Figure 1.The flowchart of the method.

Figure 2 .
Figure 2. A sample of a road profile.

Figure 1 .
Figure 1.The flowchart of the method.

Figure 1 .
Figure 1.The flowchart of the method.

Figure 2 .
Figure 2. A sample of a road profile.

Figure 2 .
Figure 2. A sample of a road profile.
2.1.1.PreprocessingCorrespondingly large data volumes and the complexity of urban street scenes increase the difficulty of creating a unified road model.Therefore, we use vehicle trajectory data (L) to section the point clouds into a set of blocks at an interval (d), as shown in Figure3.To ensure that the road in each block is as flat and straight as possible, the value of d should be set smaller under undulating and winding road conditions.

Figure 3 .
Figure 3.An illustration of sectioning road data.

Figure 3 .
Figure 3.An illustration of sectioning road data.

Figure 3 .
Figure 3.An illustration of sectioning road data.

Figure 5 .
Figure 5.A sample of elevation filtering: (a) point clouds before elevation filtering (coloured by elevation); (b) point clouds after elevation filtering (coloured by elevation); and (c) elevation histogram.

Figure 5 .
Figure 5.A sample of elevation filtering: (a) point clouds before elevation filtering (coloured by elevation); (b) point clouds after elevation filtering (coloured by elevation); and (c) elevation histogram.

Figure 10 .
Figure 10.The model of zebra crossings: (a) the distribution direction is perpendicular to the road centreline; and (b) the distribution direction is oblique to the road centreline.

Figure 11 .
Figure 11.Principal component analysis of road markings.

Figure 11 .
Figure 11.Principal component analysis of road markings.

Figure 12 .
Figure 12.An overview of the experimental data.

Figure 12 .
Figure 12.An overview of the experimental data.

16 Figure 12 .
Figure 12.An overview of the experimental data.

Figure 14 .
Figure 14.The results of road surface segmentation.

Figure 14 .
Figure 14.The results of road surface segmentation.

Table 1 .
Description of the datasets.

Table 1 .
Description of the datasets.

Table 1 .
Description of the datasets.

Table 2 .
Parameters of road marking extraction.

Table 2 .
Parameters of road marking extraction.

Table 3 .
Recognition results of zebra crossings.

Table 3 .
Recognition results of zebra crossings.