Combined Lane Mapping Using a Mobile Mapping System

: High-deﬁnition mapping of 3D lane lines has been widely needed for the highway documentation and intelligent navigation of autonomous systems. A mobile mapping system (MMS) captures both accurate 3D LiDAR point clouds and high-resolution images of lane markings at highway driving speeds, providing an abundant data source for combined lane mapping. This paper aims to map lanes with an MMS. The main contributions of this paper include the following: (1) an intensity correction method was introduced to eliminate the reﬂectivity inconsistency of road-surface LiDAR points; (2) a self-adaptive thresholding method was developed to extract lane markings from their complicated surroundings; and (3) a LiDAR-guided textural saliency analysis of MMS images was proposed to improve the robustness of lane mapping. The proposed method was tested with a dataset acquired in Wuhan, Hubei, China, which contained straight roads, curved roads, and a roundabout with various pavement markings and a complex roadside environment. The experimental results achieved a recall of 96.4%, a precision of 97.6%, and an F-score of 97.0%, demonstrating that the proposed method has strong mapping ability for various urban roads.


Introduction
High-definition mapping of 3D lane lines is widely needed for various tasks [1,2], such as highway documentation, ego-vehicle localization with decimeter-or centimeter-level accuracy [3][4][5][6][7], behavior prediction of vehicles on the road [8,9], and safe route and evacuation plan generation [1,10,11] in autonomous driving systems. Conventional maps only provide street-level information for navigation and GIS systems, while a high-definition map for autonomous driving systems requires the accurate location of every lane. 3D lane lines (lane-level maps) are one of the most important features of high-definition maps. In some GNSS-based studies [12][13][14][15][16], lane-level maps have been generated by driving along the centerline of the road (or lane) with a moving vehicle and then analyzing the GNSS trajectory. Such methods are neither efficient nor accurate because each lane must be driven in multiple times and the path tends to deviate from the real centerline. Therefore, lane mapping requires more efficient and accurate methods.
Some vision-based methods have been proposed to automatically detect and track lane markings from images or videos [17][18][19][20]. However, images and videos are susceptible to shadows cast on the road, natural light and weather conditions [21]. Because images and videos are two-dimensional recordings and lack depth information, photogrammetry-based methods based both on cars [22,23] and unmanned aerial vehicles (UAVs) [24,25] can be used to calculate the 3D coordinates of lane markings through epipolar geometry. However, successful depth measurement depends on texture. Because of the repetitive nature of lane markings, such methods are prone to mismatching and lead multithreshold segmentation followed by a morphological operation. Soilán et al. [43] selected cells with similar angles in the intensity image as a mask and applied the Otsu thresholding process over each mask. Ishikawa et al. [45] and Hůlková et al. [46] generated 2D images from 3D points to accelerate the calculation. Yu et al. [44] developed a trajectory-based multisegment thresholding method to directly extract road markings from the point cloud. In [44], interference points with similar intensities to road markings were removed by a spatial density filtering method. In contrast to the aforementioned studies, Cheng et al. [47] first performed scan-angle-rank-based intensity correction to directly attenuate the intensity inconsistency and then used the neighbor-counting filtering and region growing processes to remove interference points. For point cloud data, the combination of these two steps could improve the result of road marking extraction. However, the point cloud becomes very sparse and the quality of reflective intensity deteriorates because of the increasing range, the instrument settings, velocity of the car, and the quality of used MMS. In addition, the occlusion of other vehicles results in missing lane-marking points. In such cases, it is difficult to correctly extract lane markings using methods that depend only on point clouds. To meet the need for high-precision lane mapping, the fusion of MMS images and point clouds is necessary.
Preliminary studies on the fusion of MMS images and point clouds have been carried out. Hu et al. [48] and Xiao et al. [49] first estimated the potential road plane using the point cloud and then exploited this knowledge to identify image road features. Premebida et al. [50] generated a dense depth map using the point cloud, and combined the RGB and depth data to detect pedestrian. Schoenberg et al. [51] assigned image pixel grayscale values to the point cloud and then segmented the dense range map. Based on the findings from these studies, the fusion of MMS images and point clouds could make full use of the rich data source acquired by the MMS. To date, however, there are few studies for lane mapping. Thus, this paper aims to introduce a validation method that fuses image texture saliency with the geometrical distinction of the lane-marking point cloud to improve the robustness of lane mapping.
In this paper, a combined lane mapping method using MMS data is proposed. As shown in Figure 1, the proposed method mainly includes the following steps: 1. Road-surface extraction: After the point cloud is gridded, the sudden change in the normal vector is used to locate the curb grids on both sides of the trajectory. Road edges are fitted by the curb grids and used to segment the road points. The inconsistency of the reflective intensity of the road points is corrected to facilitate the subsequent lane extraction. 2. Lane-marking extraction: The 3D road points are mapped into a 2D image. A self-adaptive thresholding method is then developed to extract lane markings from the image.
3. Lane mapping: Lane markings in a local section are clustered and fitted into lane lines. LiDAR-guided textural saliency analysis is proposed to validate the intensity contrast around the lane lines in the MMS images. Global post-processing is finally adopted to complement the missing lane lines caused by local occlusion.
The main contributions of this paper include the following: an intensity correction method was introduced to eliminate the reflectivity inconsistency of road-surface LiDAR points; a self-adaptive thresholding method was developed to extract lane markings from complicated surroundings; and LiDAR-guided textural saliency analysis of the MMS images was proposed to improve the robustness of lane mapping.
The remainder of this paper is structured as follows. Section 2 introduces the proposed methods of lane-marking extraction and lane mapping. The experimental data and results are presented in Section 3, and the key steps are discussed in Section 4. Concluding remarks and future work are given in Section 5.

Road-Surface Extraction
As shown in Figure 2, the point cloud acquired by the MMS contains not only roads, but also buildings, trees, and road facilities such as utility poles. However, lane markings are located only on the road surface and near the vehicle trajectory. In addition, the variation of the elevation of the road surface in a local road section is gentle due to the planar structure of the pavement. To extract lane markings efficiently, this paper roughly delineates the range of the road surface based on the trajectory and elevation statistics, and then removes the points that are outside this range. For the roadside objects in the remaining points, such as flower beds, grasses, and walkways, this paper uses the normal vector to search for the road boundaries and precisely segment the road surface. The reflective intensity of the point cloud is used to distinguish the lane markings from the road background. However, the intensity of the same material in different scanning ranges is inconsistent. To eliminate this inconsistency and obtain the exact intensity contrast between lane-marking and road-surface points, intensity correction is carried out to facilitate the subsequent lane-marking extraction. Therefore, this Section takes the following three steps: (1) Preprocessing: The non-road points are quickly removed. (2) Road extraction: the road surface is extracted by the normal vector. (3) Intensity correction: The intensity of the road points is corrected.

Preprocessing
To quickly remove non-road points, a rough range of the road surface is delineated based on the dense distribution of the pavement points with respect to distance and elevation. Since the road width is limited, the pavement range on the X-Y plane can be determined using the distance from the trajectory. In general, road points in a local road section have similar elevation due to the planar structure. Therefore, the pavement elevation can be obtained based on elevation statistics.
First, the horizontal distance from the point to the trajectory on the X-Y plane is calculated, and points with distance greater than W t are removed, as shown in Figure 3. We set W t = 20 m based on the maximum width of urban roads. Second, the pavement elevation is estimated. Due to the dense distribution and similar elevation of the road points, it is assumed that the elevation histogram of the point cloud in a local road section will show a peak at the pavement elevation. To improve efficiency, a coarse-to-fine statistical method is used. The elevation histogram in meters is first plotted (Figure 4a), and the peak value is set as the coarse pavement elevation. Then, the elevation histogram in decimeters is plotted (Figure 4b), and the fine pavement elevation is obtained. Considering the possible elevation variation of the road slope, points that are 0.5 m below and 1.5 m above the pavement elevation are removed. Figure 4c,d are the point clouds before and after removing non-road points with statistical elevation, respectively. . Preprocessing based on the trajectory. The blue arrow indicates the trajectory, and the green arrow is perpendicular to the trajectory. The points outside the red line, whose distance to the trajectory is greater than W t , will be removed.

Road Extraction Based on the Normal Vector
Since an accurate road width cannot be obtained in advance, there are non-road points such as grasses and trees on both sides of the roughly extracted road points. A precise road surface should be extracted to narrow the search range for lane markings. The normal vector of the point cloud is nearly vertical on the road surface, and changes at the road boundaries due to geometric discontinuities, e.g., curbs. Therefore, this paper first locates the curbs using the sudden change in the normal vectors on both sides of the trajectory, and then segments the precise road surface using the road edges fitted from the curbs.
First, the point cloud is organized into regular grids. Then, the bounding box of the point cloud is calculated. With the bounding box, the point cloud is divided into regular grids of a certain length and width, such as 0.5 m, which is equal to the width of road shoulders. The points are then assigned to the corresponding grids based on the X and Y coordinates.
Next, we start from each grid that the trajectory passes through and search for the curb grids along the orthogonal direction of the trajectory, as shown in Figure 5. The PCA algorithm is used to obtain the normal vector of the grid, and the cosine ([0,1]) of the angle between the normal vector and the upright vertical vector is calculated. For each search, the first grid whose cosine is less than A t is considered the curb grid. The threshold A t is set empirically to 0.96 in the experiment. Finally, the road edges are fitted using the RANSAC algorithm. As shown in Figure 6a, obstacles, such as vehicles on the road, are incorrectly recognized as curbs. The RANSAC algorithm is used to eliminate the false curb grids and fit the inlier ones to road edges. 3D points located between the edges are retained as road points. Meanwhile, the obstacle curbs whose cosine are less than A t are removed. Figure 6b presents the result of road extraction.

Intensity Correction of Road Points
Due to the material difference between lane markings and asphalt or cement pavement, the reflective intensity of road markings is significantly higher than that of the surrounding road points. The intensity can be used to extract road markings. However, in addition to the object material, the intensity measured by the MMS is related to other factors, including the range, incidence angle, sensor characteristics, and atmospheric condition. Since the LiDAR data of a project are measured by the same sensor, the sensor characteristics can be regarded as constant. The atmospheric condition is negligible due to the short range of road points with respect to the vehicle. Moreover, since the height of the scanner is constant during the survey, the cosine of the incidence angle can be described as a function of the range. Therefore, the impact of the range should be eliminated.
Inspired by Kai et al. [52], a correction method is developed to eliminate the range effect by normalizing the original intensity value to a user-defined standard range. A polynomial is established to depict the relationship between range and intensity: where LI is the measured intensity of LiDAR points, R is the range, η i denotes the polynomial coefficients, N is the degree of the polynomial, and C denotes other constant factors. Coefficients η i (i = 1, 2, . . . , N) are estimated using least square regression using the road points. We set N = 3 in the experiment based on the minimum mean square error. The corrected intensity is shown as follows: where R S is a user-defined standard range that is set to the average range of the road points. Figure 7a,b are the point clouds before and after intensity correction, respectively.

Lane-Marking Extraction
Due to the planar feature of the road surface, 3D road points can be mapped to a 2D image. In the 2D image, it is efficient to apply image processing techniques to extract lane markings. Since the elevation information is necessary for the 3D lane map, a correspondence between the 3D point cloud and the 2D image is established and kept. Because of the unevenness caused by rutting and the dirt or dust coverage on the road, lane markings are mixed with noisy and interfering points. As shown in Figure 8, a low threshold value of the reflective intensity leads to considerable noise in (a), while a high threshold value leads to incomplete extraction in (b). Thus, a self-adaptive thresholding method is developed in this paper. First, the large interference region is adaptively located by a median filter, and iteratively suppressed by connected component labeling. Next, the road markings are extracted by maximum entropy thresholding. To remove false positive markings such as zebra crossings, arrows, and stop lines, shape analysis based on the geometric differences is performed. Therefore, this Section takes the following three steps: (1) 3D-2D correspondence: The enhanced road points are mapped into a 2D intensity image. (2) Self-adaptive thresholding: A binary image that contains road markings is generated.

3D-2D Correspondence
Before being mapped to a 2D image, candidate lane-marking points are selected based on the intensity statistics. First, the number of lane-marking candidate points S is estimated based on the laser scanning resolution, the length of the road section, and the width of the lane marking. Then, the road points are sorted by the intensity value in descending order, and points of the top S intensity values are kept as lane-marking candidate points.
Based on the top view of the candidate points in Figure 9, 2D intensity and elevation images are generated by dividing them horizontally into grids. One pixel of the 2D image corresponds to a grid of top-viewed points. In the experiment, according to the scanning resolution of the points, the length and width of the grids are set to 0.05 m. Since one pixel of the 2D image may correspond to multiple 3D points, the gray value of the 2D intensity image is determined by the average intensity of the 3D points in each grid. To store the elevation information, a 2D elevation image of the same size as the 2D intensity image is generated, and its gray value corresponds with the average elevation. The 3D coordinates of the pixel can be calculated by averaging the 3D coordinates of points in the corresponding grid, as shown in Figure 9. 3 × 3 mean filtering is applied to the 2D images to fill the very small holes caused by the uneven density of the point cloud. Figure 9. The correspondence between the 3D point cloud and the 2D images. On the left is a top view of the 3D point cloud, whose origin is at the lower-left corner. On the right are the generated 2D intensity and elevation images, whose origins are at the upper-left corner. The gray values of the pixel (x p , y p ) in the 2D intensity and elevation images are equal to the intensity and elevation of the 3D point (X p , Y p , Z P ), respectively.

Self-Adaptive Thresholding
For the noisy and interfering points in Figure 10a, a self-adaptive thresholding algorithm is developed. These points are iteratively located and suppressed until there is no large stretch of interference near the lane markings, as shown in Figure 10b. Finally, maximum entropy thresholding is applied to segment the road markings and a binary image is obtained, as shown in Figure 10c. First, median filtering is used to connect the noisy and interfering pixels in the 2D intensity image. Considering that the width of the lane marking accounts for approximately w l pixels and is much less than that of the interference region, the template size is empirically set to 4 · w l + 1. Second, maximum entropy thresholding and the connected component labeling algorithm are applied. A connected component with an area larger than the threshold s t is regarded as an interference region. In the experiment, we set s t = 2 · w l · l sec , where l sec is the length of the road section. Third, the regional median gray value is subtracted from the interference region. Since the intensity of the road markings is still higher than that of the surroundings, the lane markings in this region are preserved. The above steps are repeated until there is no connected component with an area larger than s t . Finally, maximum entropy thresholding is performed to segment the road markings. Morphological closing is performed to close some very small holes of the segmented lane markings, and connected components that cover less than 0.02 · w l · l l are removed, where l l is the length of the lane marking. Algorithm 1 summarizes the process of self-adaptive thresholding.  for i ← 1, n do 9: if areao f D i mask > s t then 10: In the binary image, apart from the lane markings, there remain road markings such as stop lines, zebra crossings, arrows, etc. Based on geometric features such as length, width, and orientation, the connected component labeling algorithm and morphological operations are used to remove false positive markings.

•
A single marking of a zebra crossing is similar to a lane marking, but a whole zebra crossing can be distinguished by its crossing orientation, as shown in Figure 11a,b. First, a closing operation is performed over the binary image I bin to obtain the mask image I ZC mask . Next, connected component labeling is performed on I ZC mask , and the connected components are denoted D i mask (i = 1, 2, . . . , n). Finally, the minimum area rectangle of each D i mask is calculated and the acute angle θ between the long axis and the trajectory line is obtained. D i mask with θ greater than 60 • will be regarded as a zebra crossing and removed, i.e., I bin (D i mask ) = 0. • A stop line is perpendicular to the road but connected to the lane markings, as shown in Figure 11c. Therefore, the stop line can be first separated by an opening operation and then removed by a process similar to the first step, as shown in Figure 11d,e.
where p l and p w are the tolerance of the length and width, respectively, are satisfied, D i is regarded as M j and removed (I bin (D i ) = 0). In the experiment, we set p l = 20%, p w = 30% (Figure 11f,g).

Lane Mapping
To obtain the 3D lane lines of the lane-level map, it is necessary to cluster and refine points belonging to the same lane line, e.g., points of dashed or broken lane markings on the pavement. In a local road section, this paper clusters fragmented lane-marking points based on their deviation from the vehicle's trajectory and fits them into quadratic curves. Since the point cloud becomes very sparse due to the increase in the range, and the quality of reflective intensity deteriorates at the edge of the road, this paper starts with the near-vehicle lanes and infers lanes near both sides of the road boundary step by step. A LiDAR-guided textural saliency analysis method is developed to validate the intensity contrast around the candidate lines. Finally, a global post-processing step that integrates all local extractions is used to complement missing lines caused by occlusion. Therefore, this Section follows: Considering that the lane line is not necessarily straight, the deviation is calculated along the orthogonal direction of the line segment between two adjacent trajectory points. The steps are as follows.

•
The binary image is first skeletonized to locate the center of the lane markings and reduce the computational burden.

•
As shown in Figure 12a, for each L i , a rectangular buffer B i is generated and moved along the orthogonal direction of L i with a certain step length. The step count and the number of pixels that fall within B i , which are denoted d and s, respectively, are recorded.

•
If the lane-marking pixels fall into B i , s will be larger than the threshold s t and exhibit a peak. The peak buffer is denoted P i j , where j is the number of peaks. • For all L i (i = 1, 2, . . . , n − 1), we find P i j (j = 1, 2, . . . ) with its step count d and the pixels that fall within it.

•
The DBSCAN algorithm is used to cluster all peak buffers P i j (i = 1, 2, . . . , n − 1; j = 1, 2, . . . ) according to their step counts d. The peak buffers of the same lane line will be clustered into a group due to their similar deviation from the trajectory. Figure 12b shows the clustering in an image.

Local Lane Line Inference
The point cloud becomes very sparse with the increase in the range and the quality of reflective intensity deteriorates at the edge of the road. Since several standard values of the lane width are set in the Technical Standard of Highway Engineering (2014), this paper starts with the near-vehicle lane and infers other lanes near both sides of the road boundary step by step until the curb grids are met. Based on the intensity, buffer analysis is performed around the inferred lane lines.
First, the lane width W is determined. The width of the obtained near-vehicle lane (or host lane) is calculated, and W is set to its closet standard value. Next, since the host lane is moved to infer other lanes, the curb grids located in Section 2.1.2 are selected as the boundary to delimit the movement. Finally, buffer analysis, similar to that shown in Section 2.3.1, is used to infer other lanes on both sides step by step until limit W t . As shown in Figure 13, the host lane line is moved parallel to the boundary and white pixels are searched for in a buffer that ranges from W − 0.5 m to W + 0.5 m. If white pixels exist in the buffer, the location where the white pixels match best is selected as a candidate lane line. Otherwise, the host lane line will be moved W in a parallel manner as a default candidate lane.

LiDAR-Guided Textural Saliency Analysis
In addition to point cloud data, the image data acquired by the MMS provides rich texture of lane markings and can be used to further validate the candidate lane lines. The extrinsic parameters between the laser scanner and camera can be obtained by calibrating the MMS. Based on these parameters, a potential region in the MMS image can be determined by mapping the candidate lane line extracted from the point cloud onto the MMS image. Since the lane marking is salient and linear compared with its surroundings, textural saliency analysis is applied to validate the intensity contrast in the potential region.
First, the 3D candidate lane points are mapped onto the MMS image plane. Figure 14a illustrates the relationship between the geodetic, inertial measurement unit (IMU) and camera coordinate systems. The geodetic coordinates of the 3D point [X W , Y W , Z W ] T should be mapped into the IMU coordinate system first using the geodetic coordinates of the IMU origin C W I and the rotation matrix R W I calculated by the attitude information (e.g., roll, pitch, and heading). Then, based on the collinearity equation, the camera coordinates[x, y, − f ] T of the mapped pixel can be calculated using the translation C I C and the rotation matrix R I C from the IMU origin to the camera origin, which are obtained by calibration. The formula is as follows: As Figure 14b shows, the column u and row v of the mapped pixel can be calculated by the following formulas: Considering that the width of the lane markings is approximately 0.2 m, a 1 m wide buffer around the candidate lane line points is enough for guiding the textural analysis in the image. As shown in Figure 15, because of the grayscale difference, the lane markings are the most salient in the buffer. Therefore, the existence of lane markings can be validated by textural saliency analysis. From the perspective of posterior probability, an integral histogram algorithm is used to contrast the distribution of the objects and surroundings and estimate the saliency of the pixels [53]. The candidate lane line is sampled at intervals of 0.5 m, and a group of square buffers are obtained by mapping the 1 m wide area around each sampling point to the image, as shown in Figure 16a. In each buffer, a kernel window K containing the lane marking and a border window B containing its surroundings are used to estimate the saliency of the sampling points, as shown in Figure 16b. The gray value of a pixel x in the buffer is denoted F(x). Two hypotheses are defined, i.e., H 0 : x is not salient, and H 1 : x is salient. The corresponding a priori probabilities are denoted P(H 0 ) and P(H 1 ), respectively. The saliency of x can be measured by the posterior probability P (H 1 |F(x)). Based on Bayes' theorem, we can define where p(F(x)|H 1 ) and p(F(x)|H 0 ) are the distribution of gray values in windows K and B, respectively, as shown in Figure 16b. The saliency measure reflects the contrast of the gray values between windows K and B. The width of the buffer is denoted w I , and we set the width of the kernel window w k = 0.3 · w I to cover the lane marking. Due to the slope of the lane markings in the image, the border window B is obtained by expanding the top and bottom of the window K by 25%. Since lane marking occupies a small area in the buffer, we set P(H 1 ) = 0.2 and (H 0 ) = 0.8. The window slides within the buffer to make sure that all lane-marking pixels are detected. An integral histogram is used to quickly calculate the grayscale histogram in each slide window, as shown in Figure 16c.
The saliency values in the buffer are then segmented by the threshold S t , as shown in Figure 16d. We set S t = 0.8 in the experiment. If a linear object exists in Figure 16d, the corresponding sampling point is considered valid. For each candidate lane line, we denote the number of all sampling points as n and the number of valid points as m. Considering the dashed lane markings and occlusion, the candidate lane line is considered valid if m/n ≥ 50%. An example is given in Figure 16a.

Global Post-Processing
After local lane mapping, 3D points of lane markings are sampled and validated. There are some missing lane line points due to occlusion. This paper uses the global continuity of the lane lines to complement the missing points. If the distance between two adjacent lane line points is greater than the threshold D thre , we move the corresponding trajectory points to complement the missing points as illustrated in Figure 17. Finally, B-spline curve fitting is used to smooth the lane-level map.

Experimental Data
The data in this paper were collected using a real-time 3D LiDAR scanner Velodyne HDL3 with 700,000 measurements per second, a horizontal field of view of 360 • and a typical accuracy of ±2 cm. The reflective intensity can be simultaneously measured. The navigation system integrates a GNSS satellite positioning module and an IMU. The images are obtained by a Point Grey camera with a resolution of 2448 × 2048. Therefore, the 3D point clouds, images, and trajectory data that we need in this paper can be simultaneously obtained while the vehicle moves.
The collected data are located around a bridge in Jiangxia District, in the city of Wuhan, Hubei, China. As shown in Figure 18, the data cover approximately 3.6 km of a two-way six-lane road and a roundabout with four-lane entrances. The two-way road is separated by the double yellow line or guardrail around the center and has straight and curved parts. There are diverse objects (e.g., flower beds, grasses, walkways, etc.) on the roadside and vehicles on the road. In addition, there are various road markings on the road, such as zebra crossings, stop lines, arrows, and diamond markings. Due to unevenness caused by rutting and the dirt or dust covering the road surface, the lane markings are subject to noise and interference. Moreover, the point cloud is very sparse toward the roadside, and the quality of reflective intensity deteriorates near the road boundaries.

Evaluation Method
The resulting lane-level map is compared with the ground truth that is the centerline of the lane markings and manually annotated by professional operators. The result is considered true if it falls within a 10-cm-wide buffer of the ground truth. The evaluation criteria are given as follows: where TP, FP, and FN are the lengths of the true positive, false positive, and the false negative lane lines, respectively. It is clear that higher precision means fewer false positives and higher recall means fewer false negatives.

Typical Cases
As shown in Figure 19, a typical road section is selected to illustrate the process of local lane mapping using the proposed method. The raw point cloud (Figure 19a) is a three-lane road that is separated by the guardrail and the roadside hedge, with six straight arrows, a stop line, and a zebra crossing on its surface. First, as shown in Figure 19a, starting from the trajectory (the solid purple line), the curb grids that contain the guardrail and roadside hedge are located and fitted into two edges (the dotted red line). The road points (Figure 19c) are segmented by these edges. The number of points is reduced to 40% of the raw data, which reduces the computational burden. Then, intensity correction is performed to eliminate the inconsistency caused by the range. As presented in Figure 19d, the intensity of the points at the left of the trajectory (the upper points) is significantly enhanced. Based on the laser scanning resolution, the length of the road section and the width of the lane markings, 8% of the points with higher intensity are further selected to speed up the extraction of lane markings. A 2D intensity image (Figure 19e) is generated by horizontally mapping these 3D points into regular grids, with interference (the orange part) distributed along a lane. Through self-adaptive thresholding, the interference region is iteratively located and suppressed, and a binary image that contains road markings is generated (Figure 19f). After false positive markings (e.g., straight arrows, a zebra crossing, a stop line) are removed by shape analysis, a binary image (Figure 19g) that contains only the lane markings is obtained. Finally, the lane-marking pixels are clustered based on their deviation from the trajectory, and fitted into lane lines in this road section (Figure 19h). The 3D lane lines are presented in Figure 19b. For those road sections that are separated by guardrails, only half of the road is processed. However, for those road sections that are separated by double yellow lines, the lane lines of both sides are extracted. As shown in Figure 20a, the raw data are a two-way six-lane road with roadside walkways, and the purple line is the trajectory. The points far from the trajectory have lower reflective intensity and sparse density, which makes it difficult to extract the lane markings far from the trajectory. Using the proposed method, the pavement is first extracted (Figure 20c). After intensity correction (Figure 20d), the reflective intensity of the points is enhanced. Then, a 2D intensity image is generated (Figure 20e) and segmented using self-adaptive thresholding (Figure 20f). The lane near the road boundary is missed since the points are too sparse. Therefore, the successfully extracted near-vehicle lane is used to infer this lane (the upper right lane in Figure 20g). In addition, textural saliency analysis is used to validate the inferred lane. A group of square buffers are obtained in the MMS image, and the pixels with posterior probability P(H 1 |F(x)) greater than 0.8 are considered salient, as indicated by the white pixels in the upper picture of Figure 20h. There are 120 valid sampling points of the 120 points, i.e., m/n = 100% > 50%; thus, the inferred lane is considered valid. Finally, the 3D lane lines are obtained, as shown in Figure 20b. For the double yellow lines, we first extract the line that is closer to the trajectory. By analyzing its surroundings, the other line can be obtained. For the curved road as shown in Figure 21a, the proposed method processes the point cloud piecewise based on the line segments between the adjacent trajectory points. The curb grids are first located by searching along the orthogonal direction of these trajectory line segments, as shown in Figure 21b. The pavement points are segmented by the curbs and enhanced by intensity correction (Figure 21c). An intensity image is generated using the pavement points and segmented by self-adaptive thresholding (Figure 21d). Then, along the orthogonal direction of the trajectory line segments, the lane markings' deviation from the trajectory can be calculated despite the curvature of the road (Figure 21e). Finally, the lane lines are obtained by clustering and fitting the lane-marking pixels (Figure 21f). The 3D lane lines are shown in Figure 21g.  Figure 22 shows the overall mapping results. The red lines are the inferred lane lines that are considered invalid in the textural saliency analysis. Since they can improve the continuity, they are kept for further manual check but not included in the evaluation. Three typical road types (e.g., straight, curved and roundabout) are selected and overlaid with Google Earth images [Map Data: Google, 2018 DigitalGlobe]. As shown in Figure 23, the lane-level map coincides well with the image, which demonstrates the effectiveness of the proposed method.

Evaluation and Error Analysis
As shown in Table 1, the recall, precision, and F-score of the overall lane mapping reach 96.4%, 97.6% and 97.0%, respectively. The false positives are caused mainly by the complex junction where the roundabout and the three-lane road join each other. An irregular lane with much larger width than that of the standard lanes causes false negatives. The detailed cases will be analyzed below. As shown by the point cloud ( Figure 24a) and the corresponding MMS image (Figure 24b), the junction where the vehicle enters the roundabout is rather complex, which makes it difficult to track the lanes correctly. As shown in Figure 24c, the gray lines are the incorrectly tracked lanes. In an actual project, it is necessary to first locate those junctions and annotate the lane lines manually.

Discussion of Textural Saliency Analysis
Textural saliency analysis can be used in many cases, including the long and solid lane markings near the road boundary mentioned in Section 3. There are cases when the interference is so strong and the contrast between the lane markings and the asphalt road is so weak that the method fails to extract the near-vehicle lane. If only one of the near-vehicle lane lines is successfully extracted, we infer it to obtain the other one. The width of the lane W is decided by the previous road section. As shown in Figure 27a,b, the other near-vehicle lane line is obtained and validated by textural saliency analysis. Since the slope of the near-vehicle line is larger than that of the lane lines near the road boundaries, the shape and size of the detector window are set horizontally as shown in Figure 27c. The border window B is obtained by expanding the left and right of the window K by 25%. In addition, the height of window K is set to 0.5 · w k , where w k is the width of K. After textural saliency analysis, since there are 56 valid points of the 71 sampling points and m/n = 78.9%, this line is considered valid.

Conclusions
Most current studies based on the MMS do not consider the images and only use the point cloud to extract road markings. However, the point becomes very sparse and the quality of reflective intensity deteriorates because of the increase of the range, instrument settings, velocity of the car or the quality of used MMS, which makes it difficult to successfully extract the road markings from the sparse points for those studies. The combined lane mapping method proposed in this paper can effectively solve this problem. First, intensity correction is used to eliminate the intensity inconsistency caused by the range. Then, self-adaptive thresholding of the intensity is developed to extract the lane markings against the impact of the interference. Finally, LiDAR-guided textural saliency analysis is adopted for lane mapping, especially for those lane lines that are far away and difficult to extract by methods that use only point clouds. Global post-processing complements the missing lane lines caused by the occlusion of vehicles and improves the robustness of lane mapping. The test results with the dataset from Wuhan achieved a recall of 96.4%, a precision of 97.6% and an F-score of 97.0%, demonstrating that the proposed method has strong processing ability for complex urban roads.
Future studies can use the prior map to rely less on real-time trajectories and improve the robustness of lane mapping, especially when the measuring vehicle changes lanes. Moreover, an update of the lane-level map that integrates cloud points and images is suggested for further study. A low-cost photogrammetry-based generation of lane-level maps using UAVs and cars is also suggested.