Extraction of road boundary from MLS data using laser scanner ground trajectory

2.1 Estimation of laser scanner ground tracks

This study estimated scanner ground tracks based on point density and slope in the three steps.

2.1.2 Obtaining rough scanner tracks

Points high above a road with low slope may be extracted to TY because of environmental complexity. To remove them, we segmented TY by time interval Δt, based on GPS time T. Let m be the number of segments after TY is segmented. Elevation is analyzed for each segment TY _i , where TY _i = {TY((i −1) · Δt < T _i < i · Δt)}. Those points with the most concentrated elevation distribution in TY _i are retained and marked as the i-th road feature area TYF _i , TYF _i = {TY _i (Z _P − ΔZ < Z _i < Z _P + ΔZ)}, where Z _P represents the peak value in elevation histogram of TYF _i and ∆Z (0.2 m is recommended) is the allowable elevation tolerance. Rough scanner tracks CG (red points in Figure 2) are obtained by averaging the coordinates of all points in each TYF, giving CG i = TYF i ( X , Y , Z ) ¯ .

Figure 2

The rough estimated scanner ground tracks.

In Figure 2, the rough estimated trajectory (centers of gravity of the high point-density areas horizontally distributed) swings left and right, which is not consistent with the driving rules. Therefore, a refining process is necessary.

2.1.3 Refinement of the estimated scanner ground trajectory

In Figure 2, colored strips indicate the road feature area TYF divided by time interval Δt, in which CG points (in red) swing as the point distribution varies when the vehicle travels forward. To further refine CG _i , we fitted them to an isochromatic series based on T. The CG points are first returned to the original point data based on the nearest neighbor principle (Euclidean distance). Then, the corresponding scan lines number to the CG points are detected. A test region of scan line frequency is constructed between CG _i and CG _i−1 shown in black rectangle in Figure 3.

Figure 3

Test region to obtain the scan lines CG _i belongs to.

The test region is several times wider than the average point spacing and longer than the CG point spans. The value of T in test region has slightly change in points 1–4, 5–8, and 9–11, while sharp time difference occurs at points 4–5 and 8–9. The sharp difference indicates the approximate interval time between adjacent scan lines. Let it be t _L. Then the scan line number CG _i belongs to can be derived by:

(3) SNl i = T CG i − T 1 t L + 1 2 ,

where the variable SNl _i denotes the sequence number of the scanline CG_i belongs, take the scan line of the first point of data set P as line 1.

All CG _i points from multiple measurements of laser foot points corresponding to the same scanning angle, and their time T should meet the isochromatic sequence shown in equation (4).

(4) T CG i = t 1 + ( SNl i − 1 ) ⋅ t d ,

where t ₁ and t _d represent the reliable estimate of acquisition time of the refined ground track points on the first scan line, and the elapsed time of each scan line, respectively. The point whose acquisition time closest to the refined T CG i is taken as the estimated vehicle tracks which is expressed as NPH _i .

2.2 Composition of edge blocks

The raised curb is often detected as quantitative and discrete candidate boundary points. This study combines continuous fluctuating points into edge blocks by searching forward and backward from NPH _i . Figure 4 illustrates the process of building edge blocks in four consecutive steps.

Figure 4

The process of composition of edge blocks.

As can be seen from this figure, in the first step, a moving window is created to detect upward fluctuating points (small squares in cyan in Figure 4) by height differences. Window size is expressed as the number of points contained within the window. The half width (BMW _i ) of the window can be determined by BMW i = ⌈ C h / ( J s i ⋅ sin θ ) ⌉ , where Ch and θ denote the minimum height (0.08 m) and minimum slope (30°) of the boundary to be detected, respectively. The default values of these parameters are considered according to regulations on exposure height of curb stones in Chinese city streets. In addition, window width varies with J s i ; the larger the window, the fewer points it contains. The elevation difference H diff i before and after the center point of the window is calculated as:

The point at which H diff i ≥ Ch is the upward point. We only detected boundaries in the forward direction, that is, d i + BMW i − d i − BMW i > ε , where d i + BMW i , d i − BMW i denote the distances of the last and the first point in the window from NPH _i . Tolerant ε is set to −0.1 m (to allow return of 0.1 m). The search length in Figure 4 is set to 15 m, according to conventional road design width. Two point sets A1 and A2 are formed by forward and backward searches from NPH _i in opposite directions. In the following steps, forward search from A1 is taken as an example.

The second step is to filter points by slope. We only consider the slope between central point i of the window and the next point i + 1, which is given by:

Points with Slope _i < tan(θ) are removed from the data. It is possible to have one or two continuous upward points that may not be near the road boundary. However, if continuous multiple points fluctuate upwards, they are more likely to belong to a boundary. We combined continuous upward fluctuation points into blocks in the third step. The consecutive number Bn _i must satisfy Bn i ≥ max ( ⌊ η C h / J s i ⌋ , 1 ) , where η is the correction coefficient determined by the ratio of the interquartile range to the total range of adjacent point spacing among 10 neighborhood points. The value of η depends on the rate of missing points; a value in the range of 0.7–1 is recommended.

In the last step, the elevation of NPH _i , extracted in Section 2.1.3, is used to remove blocks whose starting points are higher than the road surface.The filter is given in equation (7) and shown in Figure 4 with cyan shadows. Let (X _i , Y _i , Z _i ) represent the position of NPH _i and (BSX_j , BSY_j , BSZ_j ) denote the coordinates of the start points of the j-th block located on the NPH _i scan line. The threshold value of Z _th is taken from the greater of the two values of 0.1 m and 3% of the span according to the horizontal slope of the road (generally 1.5–2%).

As shown in Figure 4, the blue and green points are the first and rest points of blocks, respectively. After multiple constraints, there are few edge blocks on an NPH _i scan line. If no edge blocks remain, the boundary point corresponding to NPH _i is labeled as null. Otherwise, feature information of discovered blocks, including BSD i j , spacing from the starting point of the j-th block to PG _i , and BH_ij , the elevation difference in the block is extracted for the boundary tracking process. These features are calculated as:

where BEZ_j represents the elevation of the last point in the block.

2.3 Boundary tracking

Frequent lane switching is not allowed during scanning, which limits fluctuations in spacing between the road boundary and the scanner’s ground track. This relationship is reflected in a pseudo-mileage spacing map in Figure 5, in which NPH _i points are represented as red circles. Blue and white rectangles denote the starting point of each block and generated candidate edge feature points, respectively. The gray line shows the scan lines located at NPH _i . There may be more than one edge block extracted from NPH _i . However, the first block in each scan line is highly likely to belong to the roadside when the road surface is relatively smooth. Points on the map are named candidate boundary feature points, corresponding to ( ∑ 1 i L i , BSD i j ) of each boundary point, where L _i is the Euclidean distance between adjacent NPH _i . There are three processes to present the boundary tracking algorithm.

Figure 5

The pseudo-mileage spacing map.

Process 1: Extracting the initial feature segment of the road boundary

In most cases, y _i varies little in a certain range of Δx. However, there is always a segment from the map longer than 5 m that accords with a robust linear regression model of 95% confidence. These segments are found in the dataset composed of the first block of each scan line and are labeled the initial road boundary feature segments (red line in Figure 6). Inside points satisfying regression conditions are identified as boundary feature points, and end points at both sides are identified as the last identification points (red points in Figure 6).

Figure 6

Hunting zones derived by Δx and α.

In Figure 6, the initial feature segment of the road boundary is approximately in a straight line. We use this line as the direction to establish a possible area containing nearby boundary feature points at the end of the last identified feature points. The vertex angle of the possible area changes with the distance (x) from the vertex (the last identified points), whose size is determined from Process 2.

Process 2: Predicting the possible location of the next boundary feature point

The vehicle keeps moving in its original direction because of inertia. We first evaluated the separation of the vehicle from the road boundary based on the tail of recognized edge feature points (not less than 3 m) where Δy = y _max − y _min. If Δy _r ≤ y _th, where y _th is a set threshold, then the vehicle runs parallel to the road boundary and the possible location of the next boundary feature point is close to the last recognized y _i . Otherwise, the separation of vehicle from road boundary is represented as a regressed linear model. Then, y _i+1 of the next boundary feature point at x _i+1 can be predicted as:

where y = k _l x + k_b represents the regressed straight line derived from the tail of the recognized boundary feature points. Let ψ be the angle of the regressed line from the positive x-axis, then k _i = tan(ψ). The threshold y _th is set at 0.25 m, giving a maximum deviation of 5° from the road edge. It is common for vehicles to travel parallel to the road boundary. Then, the prediction y _i+1 = y _i improves tracking speed.

Process 3: Tracking boundary feature points using a hunting zone

We tracked boundary feature points using a hunting zone (black in Figure 6) controlled by a fan-shaped center angle α. As shown in Figure 6, α is opened along the predicted direction, the vertex of which refers to the last recognized boundary feature point (the last green point). We recommend an empirical model of α expressed as α = (a exp(−bΔx)) + c, where, a, b, and c are constant coefficients controlling the radian-unit fan-shaped angle. The recommended values of (0.86, 0.76, 0.12) cause α to decrease with the passing of Δ x . In practice, the search area is converted into a hunting zone parallel to the y-axis direction and is deduced as follows:

where [ y i + 1 − R , y i + 1 + R ] constitutes the hunting zone. If the hunting zone does not include a point where x = x _i+1, then the search continues for x = x _i+2. If multiple points are found in the hunting zone, the point with the smallest angle to the search direction is selected. However, as the distance between the current point and the final identified boundary feature point (the last green point) increases, the reliability of tracking deteriorates. We set a maximum search length SL_th of 15 m. Exceeding this length, the forward search process is interrupted and tracing in Processes 1–3 is re-implemented from the last identified boundary feature point. Figure 7 illustrates the segmentation of feature points in the tracking process.

Figure 7

Data segmentation in tracking process.

The red line on the left side of Figure 7 is the initial edge feature segment recognized for the first time. It divides dataset B into two segments B2 and B1, both of which contain red line segments. Starting from the end points (in red), the two sides are tracked in opposite directions. The identified edge feature points tracked in opposite directions are marked in green. If no point is marked in SL _th from the last identified point, the data are disconnected. Beginning with the last green identified point, a new search dataset B3 is constructed, on which Processes 1–3 are executed again. The initial boundary feature segment (red line on the right side of Figure 7) is found for the second time. B3 is divided into B4 and B5 segments for forward and backward search, respectively. Note that the overlapping search area B6 undergoes a bidirectional search, forward search in the first dataset B1 and reverse search in B5. This can reduce the probability of missing boundary points. Provided that the recording order of data points in B2 and B5 is reversed, the same search steps as those in B1 and B3 can be used.

2.4 Post-processing

Once the connection relations of the extracted boundary points on the scan line are determined, a rough vectorized boundary can be quickly obtained. The selection of boundary mathematical models and the extraction of all point clouds in the road are not addressed in this paper. The following sections focus on determining the connection between adjacent boundary points.

Length density Ld (m⁻¹) of two boundary points is defined as the ratio of the number of points in their spanning region to their distance. Figure 8 illustrates the span cube based on two adjacent points. In this figure, the red points are extracted edge points and the green lines connect two adjacent boundary points, cb1, cb2 (set at 0.2 and 0.1) depict the searching width in the direction of heading or departing from the road, and ch1, ch2 in Figure 8(c) indicate the height range of the cube, where ch1 is higher than the tallest point and ch2 lower than the lowest point (both set at 0.1). The connector Con describes the connection between adjacent boundary points as:

where A s and s i − i + 1 refer to the average point spacing of raw data and the Euclidean distance between the current boundary point i and the next point i + 1, respectively. φ i is the deviation angle formed by the line direction obtained by least squares linear regression from five points before and after point i . When Con = 1, the point i is connected to the next point; if not, points are not connected.

Figure 8

Span cube based on two adjacent edge points. (a)–(c) represent the search cube in axonometric, top, and side view, respectively; (d) illustrates the length of the span cube.

3.1 Data

Experimental data were collected with the SSW-IV vehicle laser modeling and measurement system developed by the Chinese Academy of Surveying and Mapping. The laser scanner is located at the rear of the vehicle on the right and the scanning plane is perpendicular to the direction of travel. Within seconds, five million measurements and 100 scan lines are collected in a typical measurement range of 1–500 m with an angular accuracy of 0.1 milliradians. Traveling at a speed of 30 km/h resulted in average spacing of 8 cm between adjacent scan lines, and 0.6 cm between two adjacent points right below the vehicle. In addition, the scan angle information is not included in the original .las file.

The datasets shown in Figure 1 include two types of road edge line: straight lines and curves. Two types of curbstones are featured, with right angle and S-shaped cross-sections. Dataset 1 (185.4 MB, 5.7 million points, 308 m long) displayed in Figure 9(a) is located in an urban area with multiple road lanes, high-rise buildings, trees, a green belt, poles, and wires. One side of the road edge is a straight line, and the other consists of broken lines with different tilt angles. Dataset 2 (279.4 MB, 4.8 million points, 300 m long) depicted in Figure 9(b) lies on a winding highway with towering mountains on one side and low-lying valleys on the other; it is surrounded by several trees, sidewalks, fences, and electric facilities. With S-shaped cross-sectional roadside stones on both sides, the edges curve with changing radii.

Figure 9

The test area rendered by elevation in (a) Dataset 1 and (b) Dataset 2.

3.2 Estimation of laser scanner ground tracks

The algorithm was performed on three sample datasets. The ground track of the scanner was estimated from two datasets at time intervals of 0.05 s, with an average distance of 0.4–0.6 m. Figure 10 illustrates the differences between the two estimated scanner ground tracks and the actual scanner position (black points). Red and blue points indicate the horizontal position of NPH before and after refinement, respectively. As shown in Figure 10, points lie near the location of the vehicle in the road, and red dots are mixed with blue. After refinement, the estimated scanner trajectory (formed by blue points) is smoother. The results show that the deviation of NPH (blue points) from the true value is not affected by road curvature, but only to the roughness of the road surface. To verify the precision and accuracy of NPH, we compared horizontal deviation from NPH to the actual scanner location. Actual location of the scanner is acquired at a frequency several times that of NPH. We extracted scanner position data with the closest data acquisition time T.

Figure 10

Comparison of estimated scanner ground tracks with the real scanner position.

The result shows that the deviation from the roughly estimated scanner ground tracks of Dataset 1 (such as Area A in Figure 10) to the actual trajectory is relatively similar in each scanning line, while the estimated tracks of the Dataset 2 (gravity centers of high-density areas) swinging around the actual ones. This phenomenon may be because of the greater roughness of the road surface of the Dataset 2, especially on the right side of Area B in Figure 10. However, the refined algorithm redistributes the estimated deviation from each scan line so that the refined scanner trajectory is very close to the actual one. In fact, the more scan lines involved in the refinement, the more stable and accurate the refined tracks are.

The estimated NPH shows strong uniformity with the actual scanner location, with a maximum value of 14.3 cm, mean deviation 2.1 cm and standard deviation 1.3 cm. For the proposed boundary tracing algorithm, the minimum width of the buffer is 40 cm, which is five times larger than the horizontal difference of the estimated points. Thus, NPH can be used for boundary tracking.

Generally speaking, the point density of road surface closer to the scanner is much higher than that of other road points. MLS systems acquire the complete points in the lane under the scanning vehicle, i.e., the core data of the trajectory estimation algorithm. These dense points of the lanes provide the basis of the correctness of the rough scanner tracks, and consequently, ensure the accuracy of the trajectory estimation algorithm.

3.3 Extraction of edge blocks

Only boundary points on the scan line of estimated tracks were extracted. Compared to extracting on all scanning lines, the number of sampling boundary points is reduced and processing efficiency is improved. Table 1 shows the integrity and the noise of the edge blocks extracted from two datasets.

The extracted edge blocks in Dataset 1 are shown in Figure 11. The red points represent the starting point of edge blocks, and the green are other points contained in the block. Left of the road, 832 ground track points (blue) correspond to 832 edge blocks (red) arranged neatly along the roadside in Figure 11(c). Locating the test area in a clean urban area with no obstructions on this side may be the main reason for this result. On the right, many vehicles are parked on the roadside, scattering and fragmenting the extracted edge blocks. As Table 1 shows, most noisy points were caused by vehicle occlusion in the road; outside the road, noise is rare. Of the 813 edge blocks extracted on the right side of the road, only seven points correspond to three edge blocks; the rest correspond to two at most. At the end of the green belt, two low edge blocks were not extracted.

Figure 11

Extracted edge blocks in Dataset 1. (a) Location of the starting point of the extracted edge blocks; (b) and (c) are enlarged maps of red and green boxes, respectively.

Figure 12 demonstrates the extraction of edge blocks in Dataset 2, which encountered the same vehicle occlusion as Dataset 1; most edge block noises were located within the road surface. Only five noise blocks were not in the road and each NPH _i corresponds to no more than two blocks in a scan line. Defined edge blocks are very reliable in finding the raised edge and have strong resistance to noise.

Figure 12

Extracted edge blocks in the winding road of Dataset 2. (a) Shows the location of the starting point of the extracted edge blocks; (b) and (c) are enlarged maps of green and red boxes, respectively.

Compared with other studies that extracted a greater number of candidate boundary points interspersed with noise [5,10,16,19], the proposed edge block detection method reduces noise as much as possible. In addition, even if more than one boundary point corresponds to NPH _i , there is not much noise on the road surface. It is very likely that the first block extracted on each scan line belongs to the road edge. Each edge block occupies a certain length, leaving multiple edge blocks on a scan further apart than the hunting zone width used in the edge tracking algorithm. As a result, there are only two cases in boundary point tracking: there either is or is not a point within the zone, which makes boundary tracking efficient. Figure 13(b) and (d) shows the position of the starting point of the edge block. The extracted boundary points (red) accurately locate the road edge. The boundary points detected in Dataset 2 shown in Figure 13(d) are not as neat as those of the urban road in Figure 13(b). This may be because the gradient of S-shaped road curbstone constantly changes and tends to be gentle near the road surface.

Figure 13

Location of extracted boundary points displayed on axonometric and side view. (a) Axonometric view in Dataset 1. (b) Side view in Dataset 1. (c) Axonometric view in Dataset 2. (d) Side view in Dataset 2.

With slope restriction, a defined boundary point has a slope greater than 30° to the next point. Alternatively, poor data quality, affected by noise and wetness of the road surface, may be responsible.

There is a common phenomenon of missing point clouds, which is reflected in sudden changes in spacing between adjacent points in Figure 13(d). On the side view shown in Figure 13(a), extracted boundary points are clustered together, whereas they are scattered in Figure 13(c). However, they remain the closest point to the boundary on this scanning line.

Except for the minimal contact area between the lower part of a wheel and the road, the outline points of vehicles have abrupt changes in elevation and distance. These changes are often geometrically backward, that is, the horizontal distance (d) from the scanner is reduced, which is not in line with the law of the defined characteristics of the boundary points. Only a very small number of edge blocks at the lower part of wheels were detected, and subsequently were ignored in boundary tracking process. Therefore, occlusions because of the parking cars on the roadside cause little trouble to the boundary points detection.

3.4 Boundary tracking

Table 2 lists the boundary points identified by the tracking algorithm. The completeness of the boundary tracking step is above 99.2%. The continuous boundary on the left side of the road in Dataset 1 was completely detected. In addition, disconnected boundaries (the left side of Dataset 2) and the branches (the right side of Datasets 1 and 2) result in the removal – during boundary tracking – of four, two, and two boundary points, respectively, as they do not conform to the main convergence direction of boundary points in the neighborhood.

The road boundary tracking results of Dataset 1 on perspective are shown in Figure 14. In the figure, blue points denote NPH, and red points the starting point of the edge blocks. Tracked points are also shown in green. Although travel direction (in Figure 14(e)) is not consistent with the boundary, the spacing between the travel direction and the boundary can be traced linearly. There is one boundary point omitted in T1 as it is beyond the search length in both directions. Three points in T2 exceed the hunting zone because of a sudden turn in the branch. This may also occur at the end of the arc green belt and road forks, but the number of such omitted points is small. Ignoring the form of the boundary, boundary points in the red circle in Figure 14(c) are extracted; these are actually on the border rather than the curb. At Area 2, depicted in Figure 14(d), the boundary is broken into very short, widely spaced parts. The quality of this point cloud acquisition is insufficient for boundary extraction. Moreover, long tracking length results in poor reliability. However, this is a rare case. The disconnection overlap area has been searched twice from two directions to detect as lots of tracked points as possible, while controlling for the risk of omission and misjudgment.

Figure 14

Results of boundary tracking in Dataset 1. (a) Candidate boundary points and (b) boundary points recognized after tracking. (c), (d), and (e) are enlargement maps of Area 1, Area 2, and Area 3, respectively.

In contrast, there is only a small amount of occlusion in Dataset 2. The vehicle trajectory is consistent with the road boundary in most cases, which guarantees reliability of the boundary tracking algorithm. Except for loss of the track at the top of the greenbelt similar to T3 in Figure 15(d), ideal boundary tracking is achieved using spacing between the vehicle and the road edge as a model parameter. The tracking algorithm is not dependent on the line shape of the road. To further verify the universality of the algorithm, a T-junction was selected for experiments. The results are shown in Figure 16.

Figure 15

Results of boundary tracking in Dataset 2. (a) Candidate edge points and (b) boundary points recognized after tracking; (c) and (d) are further enlargement maps of Area 4 and Area 5, respectively.

Figure 16

Extraction of boundary points at a T-junction. (a) The inner bend and (b) the outer bend.

Travel speed varies greatly when turning, resulting in different NPH spacing. However, the estimated scanner ground trajectory (blue points) in Figure 16 is still smooth, illustrating the universal reliability of the presented method for estimating scanner trajectory. The boundary points of the main road continue to track along the original direction (red points in the upper part of Figure 16(b)) when the vehicle turns to the branch. Boundary points on the lateral bend at the turning are quickly tracked because of aggregation of candidate boundary points. Both sides of the bend show great extraction. At the intersection of the road, where road boundaries are not continuous, the proposed method is more efficient and convenient than fitting road boundaries based on multiple models.

With short occlusion length, the correctness and integrity of boundary point extraction are guaranteed by the proximity and collinearity of boundary points on a pseudo-mileage spacing map. Rather than directly establishing a tracking model for boundary lines that may require many different models to fit different segments [21,41], the presented method uses simple linear models to track spacing between the road boundary and the vehicle. The proposed algorithm does not require the road boundary to be strictly parallel to the direction of travel of the vehicle, nor is it related to the line shape of the road boundary. Another advantage is the speed of the algorithm. The execution time is extremely short as it uses the original recording order. Comparatively speaking, calculation of the scanner’s ground trajectory takes longer than that of searching the road edge blocks and tracking the boundary points, as they only operate in scan lines located on PG. If vehicle trajectory data are available, the extraction process can be completed in real time.