Live extraction of curvilinear structures from lidar raw data

: In this paper, a general framework is proposed for live extraction of curvilinear structures such as roads or ridges from airborne LiDAR raw data, in the scope of present and past man-environment interaction studies. Unlike most approaches in literature, classiﬁed ground points are directly processed here, rather than derived products such as digital terrain models (DTM). This allows to detect possible lacks of ground points due to LiDAR signal occlusions caused by dense coniferous canopies. An efﬁcient and simple solution based on discrete geometry tools is described for supervised context in which the user just indicates where the extraction should take place. Fast response times are required to ensure a good man-system interaction. The framework performance is ﬁrst evaluated on the example of the extraction of forest roads in a mountainous area, as these objects are well marked in the DTM and hence provide some kind of ground truth. Good execution time and accuracy level are reported. Then this framework is applied to the detection of prominent curvilinear structures, which are much more diffuse objects, but of greater interest than roads in the scope of the present project. Achieved results show high potential of the proposed approach to help archaeologists and geomorphologists in ﬁnding areas of interest for future prospection using LiDAR data.


INTRODUCTION
The present work addresses the problem of the supervised extraction of curvilinear structures from LiDAR ground point cloud in the scope of interdisciplinary studies of past and present environment of a forested site located in a mountainous area.Airborne laser scanning, also called LiDAR for Light Detection And Ranging, is a 3D acquisition technique based on the emission of a laser beam swept over the measured scene and on the reception of echoes on hit obstacles.Combined with other data from on-board sensors and differential global navigation satellite systems (GNSS), the measure of the signal travelling time provides a detailed survey of the overflight landscape (Wehr, Lohr, 1999).In forest environment, the received signal is composed of multiple echoes that correspond to the successive hit obstacles, from the forest canopy, down to lower vegetation levels and, finally, to the ground itself.Therefore, a point classification is performed to separate ground points from low, medium, or high vegetation, and from other kind of specific features, such as water surfaces, buildings or wire lines.A surface is fit to ground points using optimization techniques to produce a digital terrain model (DTM) (Axelsson, 2000, Evans, Hudak, 2007), which can then be displayed using different visualization tools (Challis et al., 2011).But too sparse data may produce large approximations, that end users are not aware of (Jones et al., 2007).In particular, conifers are still a strong obstacle that impedes the laser beam reaching the soil (Popescu et al., 2002, Devereux et al., 2005), thus producing holes in the ground point distribution.In this exploratory work, classified ground points are processed rather than the DTM, in order to design and test a new approach more aware of points distribution in the raw data.Digital geometry tools are used to efficiently process the sparse point set.
This work takes place in the scope of a multidisciplinary project gathering different specialists to study past and present impacts of man activity on current landscape.Among them, archaeologists are looking for remains of man-made structures and their geographical organization, while geomorphologists aim at discriminating natural from anthropogenic reliefs.The studied terrain is a forested area rich of past human activities located in the heart of Vosges mountain in Eastern France.On-site spotting is one of the preliminary task to perform.
A LiDAR acquisition is of great help to survey forested areas, where on-site observations are often hindered by the vegetation (Sittler, 2004, Georges-Leroy et al., 2011).It provides bare soil visualizations allowing the detection of details that may be missed during on-site prospections.Moreover delineation tools are available to extract these geolocalized features and collect them into a geographic information system (GIS).More and more specialized detection tools are provided to help technical teams in this fastiduous task, using supervised or automated techniques (Sevara et al., 2016).But most often, the produced information must be confirmed or even completed by on-site validation.In the present project, a large amount of studied objects have an elongated shape: walls, ditches, ancient tracks, parcel limits, agricultural benches, etc.This work aims at providing simple solutions to extract those curvilinear structures from LiDAR raw data and to display them on DTM views.An illustration is given in Figure 1.After a presentation of available data and a brief survey of existing studies, section 2 exposes our approach and the selected methodology.Digital geometry notions used in this work are introduced in Section 3. Section 4 details the main framework designed for curvilinear structures extraction.Then an instanciation of this framework to the case of forest roads detection is presented and achieved performance on these well-marked objects is given in section 5. Section 6 focuses on the application to prominent curvilinear structures detection, and discusses possible extensions.Finally section 7 gives a conclusion and draws opened perspectives.

The LiDAR data
The Despite a mean density of 9.89 ground points/m 2 , a point repartition analysis (summarized in Figure 2) showed that the point cover is not homogeneous at all, and revealed many occluded places mostly due to the presence of coniferous canopies.This study points out the required level of robustness to data lacks, that the proposed method should reach.

Limits of existing approaches
Most of automated or semi-automated detection tools rely on the DTM (Sevara et al., 2016); this interpolated map is accurate enough for applications in well controlled environments.
In lack of efficient solutions to process the point cloud, it seems more practicable to use spatially arranged height values than sparse 3D points.In White et al. (2010), careful manual delineation of haul roads in DTM view show comparable results to accurate field surveys.Unsupervised solutions based on DTM analysis have also been developed, using a morphological classification and stochastic geometry tools to build and assemble a graph structure (Ferraz et al., 2016), or convolution neural network architectures (Salberg et al., 2017, Liu et al., 2019).But all these tools incorporate the interpolation errors.
Few approaches make direct use of LiDAR raw data.In Clode et al. (2007) None of these approaches try to directly exploit the spatial organization of available raw 3D ground points from the beginning of the process, and they strongly depend on local conditions such as surface reflectance or vegetation cover.
These works are not easy to transpose to archaeological or geological structures which are much more diffuse objects with geometrical properties that may largely vary along their layout.Some studies deal with local structures, numerous enough and geometrically distinctive to build or learn a model that can be applied to automatically detect additional instances (Toumazet et al., 2017, Trier et al., 2016).But most often, known artefacts are quite few, so that no learning basis can be used to train convolution neural networks.Moreover, surface data are generally not sufficient to surely discriminate these structures, so that on-site surveys by specialists are still necessary to collect additional information.
Recent advances of discrete geometry techniques (Klette, Rosenfeld, 2004, Couprie et al., 2019) have reached a break point, that makes them mature to address a lot of concrete problems dealing with discrete data, such as the addressed cloud of 3D points.This research field appeared fifty years ago to set the foundations of a new mathematical theory intended to the processing of numerical data produced by most of today's devices.It aims at providing computationally efficient and well controlled solutions to solve a wide range of problems.
The present work explores the use of recent discrete geometry tools to design a generic framework for extracting curvilinear structures directly from the LiDAR ground point cloud.In order to clearly bring out pros and cons of the approach, no attempt is made to filter the processed data or to exploit the interpolated DTM.A first specialization of the framework to forest roads was achieved.These objects are not really of greatest interest in the scope of the project, but as they are well marked in DTM views, they provide a good basis for the method performance evaluation.In a second stage, we set up another specialization to prominent elongated structures, in order to demonstrate the approach feasibility for detecting objects in better adequation to the project objectives.

DIGITAL GEOMETRY NOTIONS
Digital geometry notions used in this work are recalled here.
A digital straight line L(a, b, c, ν) is the set of points P (x, y) of Z 2 that satisfy : In the following, we note V (L) = (a, b) the director vector of digital line L, w(L) = ν its arithmetical width, h(L) = c its shift to origin, and p(L) = max(|a|, |b|) its period (i.e. the length of its periodic pattern).When ν = p(L), then L is the narrowest 8-connected line and is called a naive line (see Figure 3).A digital straight segment is a digital straight line restricted to A blurred segment B of assigned thickness ε is a set of points in Z 2 that all belong to a covering digital straight line L of thickness µ = ε (see Figure 4).The optimal line of the blurred segment is the covering line with minimal thickness.The thickness of the blurred segment is the thickness of its optimal line.Blurred segments can be detected in linear time by a recognition algorithm (Debled-Rennesson et al., 2006) based on an incremental growth of the convex hull of added points.At each step, the convex hull width is maintained and compared with the blurred segment assigned thickess.A directional scan is an ordered partition restricted to the grid domain G of a thick digital straight line D, called the scan strip, into scans Si, each of them being a segment of a naive line Ni, called a scan line, orthogonal to D. The directional scan is defined as the set: In this definition, the clause V (Ni) • V (D) = 0 expresses the orthogonality constraint between the scan lines Ni and the scan strip D. Then the shift of the period p(D) between successive scans guarantees that all points of the scan strip are traversed one and only one time.The scans Si are developed on each side of a start scan S0, and ordered by their distance to the start line N0 with a positive (resp.negative) sign if they are on the left (resp.right) side of N0 (see Figure 5).The directional scan is iteratively parsed from the start scan to both ends.At each iteration i, the scans Si and S−i are successively processed.
An adaptive directional scan (Even et al., 2019) is a dynamical version of the directional scan notion, that allows on-line registration to a moving search direction.Compared to static directional scans where the scan strip remains fixed to the initial line D0, here the scan strip follows a goal curve C while scan lines remain fixed (see Figure 6).An adaptive directional scan is defined by the set: where Ci is an estimate at position i of the tangent to the goal curve C. The last clause expresses the scan bounds upate at iteration i.

MAIN EXTRACTION FRAMEWORK
Based on the presented notions, a curvilinear structure extraction framework was designed, according to the following principle: LAZ points classified as "ground" are mapped into a 2D grid; a DTM-based control view is displayed to the user to let him localize potential structures; a stroke is drawn across the selected one and corresponding elements of the ground point grid are analyzed to recognize the structure section; in case of success, the structure is tracked and extended on each sides of this start position as far as possible; after post-processing validation and cleaning, the found structure is displayed over the DTM in the control view.In order to ensure user acceptability, the extraction process is required to hold in a fraction of a second, thus producing a feeling of live extraction.
Before giving thorough details on the extraction framework, we first introduce the DTM-based control view and the associated ground point grid, then the curvilinear structure model used.

The control view
It should be mentioned that the DTM-based control view is used only for two purposes: manual selection of relevant areas, and result display.Assuming that no DTM visualization technique is ideal for every application contexts and that generally several ones should be combined to guarantee a complete visual survey (Bennett et al., 2012, Štular et al., 2012, Mayoral al., 2017), we therefore implemented a simple but sufficient solution able to ensure good test conditions.It relies on a multidirectional hill-shading based on three directional light sources at 120 o from each other, one brighter at high incidence angle (60 o ), the two others at a lower angle (30 o ).The light set azimuth is let under user control to ensure a correct lighting of relevant details.
Here, accurate perception is not required from the user, and such simple visualization technique does not question the relevance of the extraction tool.Moreover, this low demanding user control is the only function ensured by the control view, so that DTM approximations have no impact in this context.

The ground point grid
The LAZ point set is organized in 500 m × 500 m tiles which cover the whole acquisition area.Classified points as "ground" are extracted and projected into a high resolution grid G that covers all the tiles and ensures fast access to the points stored in each cell.The resolution of this grid is set in accordance to the DTM resolution, that provides the input stroke P1P2.An adaptive directional scan in grid G is built from an initial scan S0 that joins the cells containing P1 and P2.The points contained in the cells of S0 are projected on line P1P2 to get an initial profile of the structure.
However this scan structure leads to an irregular point repartition along the profile (see Figure 7).Therefore, the grid resolution is set to a multiple NR of DTM resolution.The profile is then built by the set of points contained in a sequence of NR scans, thus ensuring a fairer repartition.In practice, NR is set to 5, leading to a grid cell size of 0.1 m in the present case.

The curvilinear structure model
In this work, a curvilinear structure S is considered as a surface produced by a spatially consistent extension of a characteristic noisy profile P along a 3D curve C. For example, a break-inslope profile may be seen as a succession of two noisy straight lines that form an angle α at their intersection (see Figure 8).
In that case, the noise tolerance threshold δ should be set large enough to absorb the terrain unevenness.
In order to take into account the spatial evolution of the structure S along the curve C, the profile P is declined in a sequence of adjacent profiles Pi.Each profile is extracted from the point cloud using successive scans Si.Each profile contains a list of 3D points and a geometric template, specific to each kind of structure to extract: roads, ridges, etc.The template embeds several attributes, such as height or position, in order to speed up the structure tracking in the next profile and to enforce spatial consistency.The position attribute is also used to update the scan strip ( Ci in Equation 3).At each profile position, slope and direction values of the curve C are estimated from previous profiles height and position.Null values are assumed at start.

The curvilinear structure detection framework
The curvilinear structure detection framework is depicted in Figure 9.It globally features two main stages: (i) initial detection of the structure from the input stroke P1P2, (ii) extension of the structure on each side of the initial detection.
Both sides processed ?no Then each extension step operates the following way: E1. centers the adaptive directional scan Di on the template position attribute of previous detected profile; E2. collects points of NR next scans; E3. sorts points by increasing distance to P 1 along − −− → P1P2; E4. runs a next profile tracking (specific to each structure).
The extension loop ends when the grid bound is reached or after a given threshold NF of profile detection failures if they are not due to a lack of ground points.This parameter allows to cross some obstacle, but it also increases the risk of false detections in chaotic areas.A default value of NF = 5 is assumed in this work.In areas with quite few ground points, undetected profile templates are discarded and the detection process may continue until a well covered area is met again.
Algorithms to detect first and next profiles are not necessarily identical.For next profiles, the former profile template is used to check the spatial consistency constraint, and if possible, to speed up the search of the new template.As nearby detected profiles are already validated at this stage, constraints may be released to accept some terrain roughness.

FOREST ROADS EXTRACTION
Forest roads are well marked objects in DTM views.Their delineation could provide some convenient ground truth to test the performance of the curvilinear structure extraction tool.Only main ones are already reported into geographical information systems, and this available knowledge must be updated and completed in prospection to get a more accurate position at the price of time consuming surveys on the field.Therefore, although these objects are not central in the scope of the project, their supervised extraction from the LiDAR data would yet be of great help.For all these reasons, their detection was implemented at first as a good basis for the extraction framework development and evaluation.

Forest road model
A forest road cross section, here called plateau, is modelled as a bounded horizontal blurred straight segment, with possible width varying from Dmin = 3 m up to Dmax = 10 m (see Figure 10).This model is quite similar to the one used in Ferraz et al. (2016) for forestry management purposes, and suits rather well to mountain roads.In the few flat areas found, bounds are often marked by surrounding ditches or by the hollowness due to terrain compaction.The blurred segment width is restricted to ∆H = 0.1 m, that is enough to cope with most road surface irregularities.Then, at least Nmin = 6 points are required to consider a plateau as being reliably detected.The horizontal area materializes the road plateau (see Figure 10).Because of the discrete nature of the point cloud, the bounds can not be detected exactly, but just estimated within an interval between the last point of the blurred segment and the next point.Bounds are considered as detected if this interval is less than ∆B = 0.5 m.By default, the estimated bound is this interval center.When both bounds are detected, the road width can be estimated and the profile position is set at the middle.But very often, only one bound is detected, for instance, due to the presence of an occluding tree besides the road.In that case, the profile position is defined by the detected border shifted by half the last estimated width when this information is available.
When it occurs at the first detection, new trials are made using next scans and the detection is aborted after NF failures.

First plateau detection
In initial structure detection stage of the framework (see Figure 9), the specialization of start profile detection (step D4 of Section 4.4) to forest roads proceeds the following way: 1. sorts collected points by increasing height value; 2. finds the ∆H-wide height interval containing the greatest number of points; 3. finds the longest sequence of consecutive points at this height interval in sorted points by distance to P1: it forms a horizontal blurred segment of assigned thickness ∆H; 4. checks the sequence size (≥ Nmin) and length (≥ Dmin); 5. checks if at least one of this horizontal sequence ends is correctly bounded: the distance of each end point to the next point out of the sequence should be less than ∆B.
In this algorithm, horizontal road sides at nearly the same height can be interpreted as a road plateau, thus leading to a false detection (see Figure 11).This may happen many times.Also, the proportion of those points at optimal height effectively inserted in the blurred segment is a good test to detect such situations.If the obtained blurred segment does not contain a large majority (70%) of these points, they are discarded from the height sequence and a new height interval is searched.The new blurred segment is compared to the former one and the longest one is kept.This strategy is certainly not perfect, but it succeeds in most encountered cases.

Next plateaux tracking
Optimizations based on a blurred segment detection are made for the detection of next plateaux, in a specialization of next profile tracking (step E4 of Section 4.4).Followed steps are: 1. test of available scan points count; 2. selection of scan point PC nearer to template center; 3. test of PC height wrt template height; 4. initialization of a blurred segment of assigned thickness ε = ∆H on point PC; 5. metrically balanced extension of the blurred segment on each side of PC with a tolerance of 3 successive points out of the blurred segment; 6. as soon as the blurred segment length reaches Lexp, assigned thickness pinching to the observed blurred segment thickness augmented with a security margin ∆ε, then extension continuation; 7. when the extension is stopped, test of obtained blurred segment size and length; 8. withdraw of extremal points if they belong to the bounds of the enclosing digital straight segment; 9. test of the blurred segment tilt angle ∆β; 10. if the detection failed, restart from step 3 with PC shifted 1 meter on the right, then 1 meter on the left.
Steps 6 and 8 are necessary to avoid the inclusion of points on road sides, that could possibly tilt the blurred segment.Therefore, Lexp is set to 2 m to incorporate most of the road irregularities before reaching the edges, and ∆ε to 0.1 m to tolerate additional roughness after the pinching.The maximal tilt angle is set to ∆β = 9.5 o to differentiate possibly slanted roads from local terrain trend.Constants were arbitrarily set, but small changes do not affect the detection process behavior.Steps 9 and 10 help to detect again the true road plateau after crossing a large flat area with wide plateaux or a pointless sector.

Achieved results
The set of plateaux which constitutes the forest road structure are displayed over the DTM view.Typically, hundred meters long road sections such as the one displayed in Figure 12

Performance evaluation
To evaluate the detector, a set of forest roads was selected in three distinct 1500 m × 1500 m large areas (see Figure 14 (a-b-c)): St-Mont (10 points/m 2 mean density), Fossard1 (7 points/m 2 ) and Fossard2 (6.5 points/m 2 ).In order to get a reference to compare detection results with, their central line was manually delineated on the DTM view.In accordance to on-site estimations, a fixed width of 3 m was assigned to most forest tracks, 5 m to larger mountain roads.As roads are not central objects for the project, no attempt was made to filter the detection outputs or even to estimate their width.Detected roads were drawn with successive plateaux connected together to fill the detection failure gaps, so that they could be compared to delineated roads on a pixel count basis (see Figure 13).Because the tool was not designed to detect road ends and forks, we artificially removed the plateaux beyond these limits in order to ensure a comparison on the basis of delineated roads, excluding possible derivations or extensions.All parameters were set beforehands and kept fixed during all the tests, including the profile detection failures threshold NF .Detected roads are displayed in Figure 14 (d-e-f).Figure 15 shows their detailed overlay on delineated roads in St-Mont area.
Let call L the cumulated length of selected roads, N the number of input strokes used for the detection, G (resp.D) the set of displayed pixels of delineated (resp.detected) roads, and S the count of pixels of set S. Achieved results are listed in Table 1, where the detection precision P , recall R and global performance F-measure F are given by Equation 4.

Discussion about technical performance
Time performance is about 220 ± 20 ms per 1000 m of delineated road, that suits quite well to interaction context requirements.Precision and recall results are quite satisfying if we consider the absence of any filtering.In particular, horizontal places nearby the road are largely integrated into the detection thus producing many side shifts that penalize the precision measure.In a road vectorization application, some filters could be added at post-processing stage and interpolated height from the DTM could also be used.In some sectors of the Fossard, the density is particularly low, and the structure extension often stops when no adequate profile is found after crossing long places with not enough points.This mostly affects the recall measure.
By definition, the directional scan is sensitive to large direction changes.The adaption ability only deals with side shifts.Therefore tight curves cause detection stops.If they could be detected, it would be possible to start a new directional scan.But no filter is applied and the implemented method does not handle rapid changes of the direction.More strokes thus necessary.Moreover, the performed delineation is certainly far from perfect.The fixed width assumption for forest roads is not always true, and the width estimation is just approximated.It may also affect both precision and recall measures.
The assumed road model for mountainous areas suits apparently well in most tested places, but it is no longer valid in flat places where the detection may often be attracted by a nearby terrace for instance.It is observed that some recent cutoffs are also undetected, the old laying-out being rather tracked.This drawback could probably be avoided by a verification of points classified as "low vegetation" in the raw data, but at the price of much longer execution time.
Using the proposed framework, such verifiable objects as mountain roads are efficiently extracted despite of the occurance of many areas with no ground points.In the following section, we show how it can be adapted to the case of other morphological structures more relevant to the project goals.

PROMINENT CURVILINEAR STRUCTURES
A prominent curvilinear structure is defined by the displacement of a convex profile in orthogonal direction along a 3D curve.This rather large definition covers a wide range of shapes at different scales, from small communal delimitation walls less than one meter high, up to cliffs or crests which may reach more than fifty meters high in the studied sector.In the present project, archaeologists are studying the remains of an undated decayed wall surrounding the Saint-Mont summit while geomorphologists pay a particular attention to breaks in slope which might delimit some spaces potentially arranged by men.
Our purpose here is certainly not to try to interpret these shapes, but rather to provide specialists with geometric tools aiming at easily extracting them from the LiDAR.Therefore, quite a neutral terminology, such as ridge, is intentionally used for these shapes.Here again, exploiting the spatial consistency of the profile along the curve is the key idea to detect these structures.

Ridge model
The ridge structure is defined by a set of cross sections, called bumps.The bump has a straight basis over which a volume can be determined through its center of mass and its height based on the given sequence of points.
The construction (see Figure 16) is as follows: let L be the baseline of the convexhull of the point set, A and B its bearing points, S the polyline joining all the points of the sequence, and S (resp.E) the nearest polyline vertex to A (resp.B) at distance HB from the baseline.We compute the area between the line (SE) and the polyline S, and determine an orthogonal line M to B that splits the polygon in two parts with same area.The bump center of mass belongs to M. We finally compute the bump position C at the intersection of M and S, and the bump height H as the distance between C and the line (SE).Note that point C does not necessarily coincide with the bump top T , defined as the highest polyline vertex to the baseline.A volume property such as the center of mass is much more a stable measure than the top, a surface feature submitted to large shift caused by terrain irregularities.C is thus better suited than T to ensure a stabilized tracking of the ridge structure.Initially intended for ridge characterizations, the base height parameter HB is not yet used, and is set to 10% of the bump maximal height, just to cope with global terrain roughness.

Bumps detection and tracking
The size of the ridge structure to detect is determined by the length of input stroke P1P2.The bump detection algorithm realizes the previous construction from the sequence of sorted points by distance to P1.It produces a template that contains the point C, the bump height H and the base points S and E. The baseline construction is optimized using Algorithm 1, where the function pointsOver(S,L) builds the set of points over the line L in set S. The other operations of the detection method are trivial.
Algorithm 1: getBaseLine: provides base line bearing points.input : A sequence P of sorted points by horizontal distance with at least one point over the line between first and last points.output: (D1, D2) → Bearing points of the baseline B ← line (first (P), last (P)); S ← highestPoint (P); P1 ← subSequence (P, first (P), S); P2 ← subSequence (P, last (P), S); D1 ← lowestPoint (P1); Next bump tracking realizes the same operations than the first detection.No particular optimization based on previous detection was found with this very simple and multi-purpose model.Therefore the template is only used for spatial consistency checking through the use of rather large tolerance values: the tolerance to point C shifts is set to ∆C = 0.5 m, and the tolerance to height variations to ∆H = 0.25 m.

Achieved results
As no real ground truth is yet available to evaluate the ridge detection, mostly visual tests were performed on well marked salient features in the DTM.The analysis of obtained bump profiles indicates a good consistency between successive scans.Moreover, the structure sinuosity is correctly surveyed, as long as direction changes are not too strong (around 45 o limit).But of course, this is not more than a visual subjective appreciation.
To corroborate these observations, we performed detection of the edge of Saint-Mont summit platform, which position is well identified.Because it features sharp angles, four input strokes were necessary to extract the whole circumference, but a correct delineation could quickly be obtained (see Figure 17).

Discussion about suitability to project needs
These results show good potential of the proposed framework to help specialists in finding areas of interest for future prospection using LiDAR data.Beyond simple structure localization, some volume quantification may also be delivered for ridge shapes whose baseline matches the local terrain trend.Notice also that the proposed framework could be easily extended to the detection of hollow structures by applying a symmetry around the baseline to the ridge model.However no such structure has been surveyed by archaeologists in the studied area untill now.
The last test shows the present limits of the ridge model used.In this case, a more complete specification of the wall profile may help to better discriminate the structure from nearby ridges.All these results must now be confirmed through detailed evaluations on concrete tasks in close cooperation with the archaeologists and geomorphologists involved in the project.

CONCLUSION AND FUTURE WORKS
This paper introduced a new framework based on recent discrete geometry tools for live extraction of curvilinear structures from raw LiDAR data as an alternative to standard DTM processing.In particular, ground point density disparities are taken into account with the new approach.To evaluate the framework performance, it was first instanciated for the detection of forest roads, whose position can easily be identified on DTM views.In regards of the simple road model used and in absence of any filtering of the extracted data, precision and recall performance achieved are quite satisfying.A long section may be extracted from a single user input and the response is immediate so that good interaction is ensured to fulfill given task requirements.A remaining challenge is the detection of tight changes of direction, that could help to reduce the amount of necessary user inputs.
This exploratory work was realized in the scope of an interdisciplinary project on past and present environment studies.Therefore, the framework was also applied to the detection of curvilinear prominent structures, and tested on the edges of a large flat area and on parts of decayed wall remains.Visually consistent detections are obtained, inciting us to develop more thoroughly this approach for effective needs related to the project goals.
Future works consist in experimental validations where detected structures will be compared to on-site surveys for finer validations in close cooperation with field specialist teams.In parallel, we intend to design a fully unsupervised detection tool based on this approach.For this purpose, visual indices could be extracted from the DTM view to automatically select relevant input strokes for the detection algorithm.Another development deals with the detection of tight direction changes using a finer estimation of the structure orientation. 50m

Figure 1 .
Figure 1.Extraction of a huge ancient stone structure (on left) from a section in the cloud of ground points (on right).

Figure 2 .
Figure 2. Ground points repartition: percentage of 1 m 2 cells with at most N ground points.

Figure 5 .
Figure 5.A directional scan: start scan S0 in red, odd scans in green, even scans in blue, bounds of scan lines Ni with dashed lines and bounds of scan strip D with bold dashed lines.

Figure 6 .
Figure 6.An adaptive directional scan.The scan strip is dynamically fit to the curve position, so that scans are shifted to continuously cover the goal curve C.

Figure 7 .
Figure 7. Profile extraction in a point grid with DTM resolution (top) and in an oversampled point grid (bottom).

Figure 9 .
Figure 9.The main detection framework.The initial detection performs the following steps: D1. creates an adaptive directional scan with S0 = P1P2; D2. collects points in NR first scans; D3. sorts points by increasing distance to P 1 along − −− → P1P2; D4. runs a start profile detection (specific to each structure).

Figure 10 .
Figure 10.Forest road section model: the road plateau lies between the displayed blurred segment and next points highlighted by the vertical strokes.

Figure 11 .
Figure 11.Risk of false detection due to horizontal road sides at the same height: the dashed height interval may contain more points than the plateau but it is split in two disconnected sets.

Figure 12 .
Figure 12.Example of road extraction displayed in white color; (a) input stroke in black; (b) section with no ground point.

Figure 17 .Figure 18 .
Figure 17.Ridge detection: edge of Saint-Mont platform.We also tested the detection of a more diffuse structure: decayed remains of an ancient wall that forms a long detectable relief ranging up to 1 m high.Parts of this structure were surveyed by archaeologists and the achieved detection conforms well to the wall laying-out.The lower part of this wall (see Figure 18 bottom left) was rather easily extracted as it stands out distinctly in this area.Only one stroke was used here.The upper part extraction was more tedious because the structure lies few meters back from a terrace edge.Unfortunately much less ground points are available in this area.The wall extraction could only be ensured by some trials before finding a correct length for the input stroke to avoid the extraction being captured by the sharper terrace border.
LiDAR acquisition campaign took place in December 2018.It covered a large part of the Fossard mountain, in an area delimited by Remiremont, Eloyes and Bémont (Vosges, France), A Titan DW device was used, with 300 kHz pulse frequency and 37 Hz scan frequency.Flight speed was 58 m/s, at 1150 m height for a slope ranging from 400 m in Moselle valley up to 800 m at Fossard top.The scan angle was set to ±21 o to ensure a 50% swath overlap.Obtained data are composed of 162 tiles with 500 m × 500 m size: 92 tiles (about 23 km 2 ) cover the mountainous sector, and the rest corresponds to surrounding urbanized plains.A set of classified points stored in LAZ format (compressed version of ASPRS laser file format) and a 1000 × 1000 pixels resolution DTM are provided with each tile.

Table 1 .
Forest road performance evaluation results.