High-Quality Contour Line Generation from LiDAR Point Clouds for the Area of Forests

A methodology for both accurate and smooth contour line generation from Light Detection 1 and Ranging (LiDAR) point clouds is proposed in this paper. In order to improve the accuracy of 2 contour lines in the area of forests, constrained triangulation networks with break lines are then 3 constructed to generate contour lines. In break line extraction, a bi-threshold method for edge line 4 detection is used to extract both complete and reliable break lines. A point clouds elevation adjustment 5 with constrain of break lines and an interpolator considering a contour interval is proposed to improve 6 the smoothness of contour lines. The proposed interpotator is also can avoid contour line intersection 7 when contour lines are interpolated. Statistical parameters and shape index are then used to evaluate 8 quantitatively the accuracy and smoothness of the resultant contour lines, which fill in the blank 9 of contour lines evaluation in theory. The experiments show that high-quality contours in terms of 10 smoothness and accuracy can be generated from LiDAR point clouds. 11


Introduction
Traditionally, contour lines can be generated by photogrammetric methods [1][2][3].However, despite their common use, these methods are quite expensive and low efficiency.Furthermore, the accuracy of the contour lines they generate is not high enough in some areas such as densely vegetation covered forest, deserts, and coastal zones, among others.In the past decade, a new generation of active surveying technique, airborne light detection and ranging (LiDAR), has been utilized due to its high accuracy in 3D spatial data acquisition, high efficiency, and relatively lower cost.One of the primary use of the LiDAR system is for DEM acquisition with a reasonably high accuracy [4].However, a critical step in contour line generation from LiDAR data is data filtering, which involves distinguishing ground from non-ground points.Almost all filtering algorithms only depend on changes in slope and elevation [5][6][7][8] , because LiDAR data are recorded in terms of 3D geographic coordinate values rather than the spectral information of ground objects.Although most modern LiDAR systems record the intensity value of a returned signal, these values are un-calibrated.More importantly, distinguishing between ground and non-ground points only by intensity values is difficult [9,10].Therefore, the ground points in the vicinity of break lines are usually classified as non-ground points.This is due to the sudden changes in elevation or slope over these areas [5,11,12], which leads to the distortion of the contour lines around them.In addition, if contour lines are extracted directly according to the point clouds that have been filtered, the resultant contour lines are prone to be saw toothed, which should be avoided because this destroys the smoothness of a contour line.Therefore, for a smoother contour line generation, the interpolation of the contour points extracted from the point clouds is necessary.This paper presents a new approach for smooth contour line generation from LiDAR point clouds with a reasonably high accuracy.The main idea behind the approach is the consideration of break lines as the constraints to overcome the distortions of the resulting contour lines.After the break lines are extracted, they are used as the constrained vectors to establish constrained irregular triangular networks [13,14].
Meanwhile, point clouds are elevation-adjusted with consideration of the break lines.Smooth and high accuracy contour lines can then be generated by this method.The next section presents a method for extracting break lines.Then the method of contour line generation with the constraint of break lines, with which smooth contour lines can be generated, is proposed.The method to evaluate the quality of contour lines is proposed in Section 4. This is followed by experiments showing the different methods for break lines extraction and the resultant contour lines.The discussion and the conclusions are described in the final section.

Extraction of break lines
Break lines are usually located along the steep zones, so this makes them easily detected by commonly used edge detectors.However, traditional edge detectors usually adopt a single threshold to distinguish edge pixels from non-edge ones, such that the extracted edge lines are either broken into too many small fragments, or too many pseudo-edge lines are formed if the threshold value is not accurately determined.Considering that a high threshold value can extract more reliable edge lines, whereas a smaller threshold value can extract more complete ones, a bi-threshold value method for extracting both reliable and complete edge lines from point clouds is proposed.The complete workflow for the extraction of break lines is shown in Figure 1.As discussed in the previous section, break lines are always located in zones where the elevation value changes dramatically [15].Therefore, borrowing edge detectors from the image processing field for break line extraction from the DEM data set is reasonable.In this paper, the LoG edge detector [16] as shown in Eq.1 is used to extract the approximate break lines.LoG is a type of second-order gradient edge detector which combines the Gaussian smooth filter and the Laplacian sharp filter into one filter where c and r refer to the number of column and row, respectively, and s = σ Gaussion .When LoG is performed for edge detection, the initial step is to obtain those pixels with zero Digital Number (DN) values after the grid DEM is convoluted by LoG, which is closely related to the threshold of gradient variance.A small threshold leads to many false but complete-edge lines, whereas a large threshold results in segmented and incomplete but less-false-edge lines.A bi-threshold strategy is adopted to combine the merits of both small and large threshold values, and the approximate edge lines are extracted by this strategy.The details are as following.

A small threshold value is used to extract complete break lines
A second-order gradient image can be obtained after the LoG operator convolutes the grid DEM.
Edge pixels correspond to those with zero DN values.Therefore, closed and connected schematic edges can be extracted by delineating pixels with zero DN values.However, large connected areas with zero DN values can be delineated in the gradient image if the point clouds are acquired over flat or very small-slope areas.This is because the gradient tends to be zero in these areas.Although complete edges can be extracted, many of these are pseudo ones because they are actually pixels where very small gradient changes occurred.A true break line should be an elongated edge.Therefore, the eight-neighboring connectivity algorithm is used to cut off those edge lines with small connected domains.In this manner, schematic break lines are successfully extracted.

A large threshold value is applied to extract high-confident edges
Large variations in gradient occur along the vicinity of a true edge, so a large threshold value is applied to extract high-confident edges.However, small variations in gradient normally exist around some true edges.In this case, a large threshold value can lead to fragmented break lines.

Combination of break lines by small and large threshold values
A small threshold value extracts more complete break lines, whereas a large threshold value detects more confident ones.Combining them results in the extraction of complete and confident break lines.Let l b be the edge pixels detected by a large threshold and l a be the edge pixels detected by a small one.They are superimposed.If the overlapping ratio in Eq.2 is greater than the predefined value, then l a is identified as a break line.

Contour lines generation with the constraints of break lines
To generate high accuracy and smoothness contour lines is always purchased by cartography community.Three strategies are considered in order to improve the accuracy and smoothness when contour lines are generated from LiDAR point clouds: (1) high accurate break lines are used as constraints for contour line generation; (2) point clouds are filtered with the constraint of break lines; and (3) contour lines are smoothed with the constraint of contour interval.

Break lines constrained triangulation networks
In the paper, break lines are used as constrained lines to be interpolated in the triangulation networks constructed from filtered point clouds in order of expressing more micro-geomorphologic details.To do so, mainly two algorithms can be employed, that is, Sloan's [17] and Floriani's [18].In Sloan's algorithm, constrained and unconstrained points are not distinguished.Constrained lines are inserted in the adjustment process followed.Floriani's algorithm is more flexible compared to Sloan's.It firstly constructs triangulation networks only by unconstrained points, and then constrained lines are inserted into the initial networks.In Sloan's algorithm, a constrained line is inserted by interchanging diagonal lines continuously, which requires the quadrangle with interchangeable diagonal lines to be strictly convex.This requirement is hard to satisfy in reality.While the core of Floriani's algorithm is Delaunay triangulation on simple polygons, the implementation of the algorithm is much simpler and straightforward.Therefore, Floriani's algorithm is adopted in the paper.The scheme indicating Floriani's algorithm is shown in Figure 2, while Figure 3 shows a resultant contour line with break lines constrained.
break lines

Point clouds elevation adjusting with the constraint of break lines
Since constrained triangulation networks constructed contains each break line(see Subsection 3.1), topographic structures can be preserved well.However, if such a triangulation networks is used for contour lines tracing directly, some contour lines would be fragmented or even saw-toothed.This is mainly caused by the following two factors: (a) the ground points obtained from filtering process contain noises, that is, non-ground points, which could generate fragmented contour lines.(b) though high density point clouds are necessary for expressing details of micro-geomorphology, it also causes saw-toothed effect of some of the result contour lines traced from such a dataset.A strategy to overcome these shortcomings is to adjust the elevation values of the ground point dataset in order to make it more suitable for fitting more smooth ground surfaces, while maintaining reasonably high accuracy.Especially in the vicinity of break lines, the result contour should be smoothly connected while preserving geomorphologic properties as much as possible.Bearing all these in mind, a break lines constrained point clouds adjusting method is proposed in the paper, whose details are described in the following.
Step 1: set accumulation error s = 0, for each point which is not located in a break lines, do the following: suppose the current point to be adjusted is P(as shown in Figure 4).Search all the neighboring points A, B, C, D, E, F within a given range, and make them as reference points.If a break line exists between P and a reference point, then find the intersection point between the break lines and the line connecting P. If they are intersected, eg.PA intersect break lines on Point A , then Point A do not set as reference points.That 's because A may be located below a cliff while P is located above the cliff of A above while P below.Therefore in this case, the reference point A can be involved in the process of adjustment.Moving surface fitting algorithm is used for adjustment in the paper.
Now suppose there are m reference points around P. Then for each point, we have the error equation: , where xi = p i (x) − P(x) and ȳi = p i (y) − P(y).For m points we have the error equations in matrix format shown in Eq.4.
where After standard adjustment computation, coefficients vector can be expressed as following: from which the elevation value of P is calculated as Figure 4.
Calculate the deviation ∆h between the estimated elevation and the elevation of check points.
The new elevation of a check point is re-calculated by Z p = (Z p + Z p )/2 and use the new value for checking if further check is required, record accumulation error s = s + ∆h.
Step 2: If s < c, then the fitted ground surface is smooth enough, stop the iteration.Otherwise, go to Step 1 for the next iteration.c = n * r, where n is the number of points, r is a predefined constant.A smaller r results more smooth ground surface.

Contour lines smoothing with the constraint of contour interval
Contour points can be extracted from the point clouds dataset after it is elevation adjusted.
However, interpolation between contour points is a must in order to generate smooth and continuous contour lines.Many interpolators have been used in literature, among which commonly used are: linear interpolation, piece-wise cubic polynomial interpolation [19], tension spline function, cubic B spline interpolation, etc.Though all these interpolators can smooth the result contour lines in some sense, none of them performs perfectly.For instance, self-intersection or inter-intersection of contour lines cannot be avoided completely.A novel tension spline function interpolator is proposed in order to guarantee smooth but no intersection contour lines, which the basic idea is that tension coefficients are adaptively setup in according to contour intervals.
Traditional method for determining whether neighbouring contour lines are intersected is calculation the topological relationship after all the contour points are connected with each other.Since there are large number of contour points in smoothed contour lines, this manner is computational expensive.We proposed a simple but effective method for the determination of whether two neighbouring contour lines are intersected.Combined with the coefficients adaptive tension spline interpolator, the following is the description of the details of the algorithm: coordinates of the points and estimate their elevation values H = h 1 , h 2 , ..., h n by linear interpolator through inserting them into the triangulation networks constructed for extracting contour points.The deviation between H and h is denoted by ∆H = ∆h 1 , ∆h 2 , ..., ∆h n , where ∆h i = h − h i .In order to guarantee the contour not intersect with its neighbors, we first give a deviation threshold ∆ < d/2, and then check the deviation point by point.Any point with ∆h i > ∆ is viewed as the contour having the trend to intersect with its neighbors, since its curvature is too large.Such a contour line must be smoothed by using larger tension coefficient σ = σ + ∆σ and repeat the procedure until all the points in the smoothed contour line satisfy ∆h i > ∆, or if the contour lies in a too steep zone to satisfy for ∆h i > ∆, then σ > Maxσ and stop the iteration instead, where Maxσ is a predefined constant.For each contour line, the workflow of the algorithm is described in the following.
Step2: calculate all the planar coordinates of points in the smoothed contour line.
Step3: estimate the elevation values of the points by inserting them into the triangulation networks constructed by structural lines constrained filtering point clouds.if ∆h i > ∆(∆h i ∈ ∆H) and σ ≤ Maxσ, then increase σ = σ + ∆σ , go to Step 2, otherwise, the current σ is used as the tension coefficient of the contour line.

Evaluation of the quality of contour lines
Quality evaluation is an indispensable step in the generation of contour lines from the LiDAR data set.The traditional method for quality evaluation of contour lines is based on ground truth, a method similar to other spatial data quality assessment approaches which is labour consuming in terms of collecting ground truth.Furthermore, the method can assess the quality only by point-based truth.Collecting ground truth corresponding to each point of a contour line is impossible.In this paper, statistical and shape-index based models are proposed for quality evaluation.The former is used for accuracy assessment, whereas the latter is used for smoothness measurement.

Accuracy assessment of contour lines by the statistical-based method
After contour lines are generated by the process described in this paper according to the given map scale, the following steps are proposed to evaluate the accuracy of the result: (1) Each contour line with a 2 mm interval is sampled, in which the interval is the length measured in the map sheet.The real coordinates of the point being sampled should be converted from the interval and the given map scale.The sampled points is set as Point Set I.
(2) A triangulation network is constructed with Point Set I.
(3) The filtered LiDAR data set is input, which can be taken as ground check points, and the elevations of the check points are estimated by linear interpolation using the triangulation network constructed in Section 3. 3 Step 2. The elevation deviations between the interpolated points and the check points are then calculated.
(4) The following statistical parameters are calculated: maximum, minimum, average error, and root mean square error (RMSE).

Smoothness measurement of contour lines by the shape-index based method
The concept involved in smoothness measurement based on shape index is that the smoothness of a curve increases with its chord length, whereas it decreases with its curvature.Given contour points B, A and C (shown in Figure 6) and a triangle constructed by them, the larger the A and the longer the chord length a , the smoother the curve.From a local point of view, A alone is enough to measure local smoothness.However, chord length a effects the overall smoothness of the curve.From this point of view, the formulafor smoothness measurement is proposed as Eq.6.For a whole contour line which is segmented by contour points, the ratio of the summation of smoothness degree by each segment to the summation of each chord length is used to evaluate its smoothness measurement as Eq.7.
where m denotes the number of triangles constructed by the contour points of the line.If the contour line is a closed one, m is equal to the number of contour points, whereas m is equal to the number of contour points minus 2 if it is unclosed.

Experiments and analysis
The procedure for contour line generation from the LiDAR data set as described in this paper has been validated by a real data set collected by Leica ALS50, the second-generation LiDAR system.The data set covers Changyang area, Hubei Province, China, a hilly area with a relative relief of more than 1000 m.The main system parameters for data collection are as follows: minimum flying height of 1300 m above average ground, scanning angle of 49 degrees, and an overlapping rate of 30%.Spacing along the flight is around 1.3 m, whereas it is 1 m across the flight.The results are in 1.2 points per square meter point density of the final point clouds.The total data volume is around 2.3 gigabytes.

Experiment 1. Break lines extracted using the bi-threshold value method
The experiment tests the efficiency of extracting break lines using the bi-threshold method.Point clouds are rendered according to their elevation values, as shown in Figure 7a.The area covers several terraced fields.
The procedure for this experiment is as follows.
(1) Filter by progressive TIN densification and then generate gird DEM by linear interpolation.
(2) Convolute the grid DEM by LoG operator, with the size of the convolution set to 25 × 25 and the Gaussian variance set to 2.5.
(3) Extract edge lines using a small threshold value (5% of the average gradient).The result is shown in Figure 7b.Then those edge lines with small connected domains are cut off, as shown in Figure 7c.
(4) Extract edge lines using a large threshold value (75% of the average gradient).The result is shown in Figure 7d.
(5) Extract schematic break lines by combining small and large threshold values.In the experiment, the overlap ratio(see Eq. 2) is 0.4, and the schematic break lines are shown in Figure 7e.From Table 4, smoothness evidently increases with a corresponding decrease in the map scale.
In a small map scale scenario, the difference in smoothness between with adjustment and without adjustment becomes negligible.

Experiment 4. Relationship of accuracy and smoothness
In areas with dramatic topographic changes, the accuracy decreases with a corresponding increase in smoothness if the contour lines are generated for a large map scale.An experiment is then designed to test the relationship of accuracy and smoothness.A terrace field is selected to represent an area with a sharp topographic change.Table 5 shows that a larger r signifies little adjustment in the point clouds.Although the resultant contours can better approximate the real topography in this case, the smoothness of the contour becomes low.

Discussion and conclusion
In the paper, a methodology for both accurate and smooth contour line generation from the LiDAR data set is described.First, highly accurate break lines are extracted.Second, filtered point clouds are then adjusted by the constraints of these break lines.Finally, constrained triangulation networks are constructed.In break line extraction, a bi-threshold method for edge line detection is used to extract both complete and reliable break lines.In contour line interpolation, an interpolator, with consideration of the contour interval, is proposed to overcome contour line intersection.Finally, statistical parameters and shape index are adopted to evaluate quantitatively the accuracy and smoothness of the resultant contour lines.The experiments show that in areas with serious topographic changes, the contour lines generated from LiDAR points after adjustment with the constraints of structural lines are more accurate than those from LiDAR points after adjustment without the constraints of structural lines.
In our proposed break lines extraction strategy, if the size of a grid mesh is too small, i.e., if the resolution of the grid DEM is too high, then more time is consumed for the extraction of break lines than low resolution.This is mainly due to a piece-wise construction algorithm is used to refine break lines extraction.This is an iterative process and is thus a time-consuming task, especially when a small step is adopted for the iteration.To overcome this shortcoming, the following are proposed in future

Figure 1 .
Figure 1.The process of extracting break lines.

Figure 2 .Figure 3 .
Figure 2. The procedure of constructing triangulation networks constrained with break lines.(a) Ground points and structural lines, (b) Construct D-TIN with the key points of break lines, (c) Insert break lines, and (d) Construct D-TIN constrained with break lines.

Figure 4 .
Figure 4. Select reference points for adjustment.Point C, D, E, F can be set as reference point for processing adjustment, while A, B can not.

Figure 5 .
Figure 5. Contour lines smoothing with the constraint of contour interval.

)Figure 6 .
Figure 6.smoothness measurement of contour line by shape-index based method.

Figure 9 .
Figure 9. Relationship of accuracy and smoothness.

Figure 9
Figure 9 shows the relationship of smoothness and accuracy of the contour lines.Evidently, RMSE increases with a corresponding increase in smoothness.Further, by comparing the solid line with the dotted line, the contours from cloud points adjusted with the constraint of structural lines can maintain more topographic features compared with those from cloud points adjusted without the constraint of structural lines.

research undertakings: ( 1 )
identification of the empirical optimal parameters of the size of the grid mesh and iteration step in terms of the efficiency of the program, and (2) further development of the programming method from the traditional single core CPU to the multicore CPU environment, or the use of another parallel programming environment such as cluster or grid computing to improve the running efficiency of the program.Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 1 February 2018 doi:10.20944/preprints201802.0003.v1

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 1 February 2018 doi:10.20944/preprints201802.0003.v1 that
smoothens the data set first and then detects the edge lines.It reduces the influence of noise for edge detection.log

Table 3 .
RMSEs of contour lines in different scales.

Table 3 ,
Method I: Contour lines are generated from points without elevation adjustment.Method II: Contour lines are generated from points with elevation adjustment but without the constraint of structural lines.Method III: Contour lines are generated from points with elevation adjustment and with the constraint of structural lines.

Table 4 .
Smoothness in different scale modes.
The smooth parameter c is tuned to generate different contour lines with different smoothness and accuracy values, where c = n × r (see details in Section 3.2).Contours with different smoothness values are generated by adjusting parameter r, and then the RMSE and average smoothness values are calculated.The results are shown in Table 5.The variation in the trend of RMSE and average smoothness (AS) is shown in Figure .

Table 5 .
Accuracy and smoothness in different processes.