Module LineModeler
See also
opals::ILineModeler

Aim of module

Models 3D structure lines based on a point cloud (ODM) and 2D line approximations.

General description

The 3D point cloud (inFile) in ODM file format and 2D approximations of the structure lines (approxFile, valid formats: ESRI shapefile, (Binary) Winput) constitute the main input for Module LineModeler. The resulting 3D structure/break lines are stored in a vector data file (outFile) either in ODM, ESRI shapefile, or (Binary) Winput file format (oFormat). By default, the entire 3D point cloud is used for structure line modeling. However, the input points can also be restricted using an appropriate filter string. For processing only a subset of lines a list of line ids or a single textfile containting the respective line ids can be specified (processId). Please note, that for line approximations provided in Winput file format, the is correspond to the structure number, whereas else (i.e. Shape and ODM) the consecutive line numbers are used. The same strategy can be applied to skip processing of certain lines via parameter ignoreId.

The general 3D modeling strategy of Module LineModeler is based on the intersection of local geometry primitives (plane, cone, or polynomial cylinder). The length and width of the surface elements to the left and the right of the 2D approximation can be set by the parameters patchLength and patchWidth. Please note, that a minimum and a maximum length can be specified denoting the smallest and largest longitudinal extension of the local surface elements. For the patch width either one value (symmetrical width to the left and the right of the 2D line approximation) or two values (1st=width left, 2nd=width right) are expected. Although usually a constant patch width is employed for modelling the entire line network, the maximum width is, however, locally constrained by the neighbouring 2D line approximations to ensure that a local surface element covers a smooth surface (see Figure 1). The patch length is automatically adjusted by Module LineModeler based on the curvature of the 2D approximations (i.e. the straighter the line progression, the longer the patch). In addition to the patch length and width, lower and upper bounds for the number of points within a single geometry primitive can be set (pointCount).

Fig. 1: Data selection strategy of consecutive patches

Once the points to the left and the right have been selected, local geometry primitive fitting is performed. It is important to know that both geometry primitives (left and right) are estimated in an integrated least-squares adjustment minimizing the vertical point-to-surface-patch deviations while considering the line approximations as additional observations. The latter helps to stabilize the results in awkward point configurations. Currently, three different patch pair fitting methods are implemented. The standard model is the plane pair fit as shown in Figure 2a. The actual result of each patch adjustment is a representative point P, the tangent t, the intersection angle, and several accuracy (σ) values. To avoid extrapolations the representative point P is computed by intersecting the two patch planes with a vertical plane that contains the center of gravity (COG) of the left (SL) and the right (SR) point set.

Fig. 2a: Standard plane pair fit

In case the 2D line approximation is strongly curved within the current patch, the modeling strategy computes a plane pair fit and a plane-cone fit. As can be seen in Figure 2b the used cone model is a circular cone with a vertical axis. The results of the plane-cone fit are selected if the σ0 a-posteriori is lower than the σ0 a-posteriori of plane pair fit. To compute the representative point P a vertical plane containing the COGs of both point sets is computed and intersected with the adjusted plane and cone. In general, there are two solutions for P which is why the distance to the approximation line is used to select the correct one.

Fig. 2b: plane-cone fit

In cases were the surface is curved perpendicular to structure line but straight in line direction, plane and cone primitives sometimes deliver unsatisfactory results. For such situation a 2nd order polynomial cylinder model was developed as shown in Figure 2c. From a geometrical point of view, the model consists of two cylindrical surfaces that are formed by the motion of two parabolic curves (the generators) across the intersection line (the directrix). Testing the statistical significance of the quadratic terms helps reducing over parameterisation. The quadratic term is deactivated if not statitically significant, and the adjustment is repeated in this case. Hence, if both quadratic terms are deactivated, the model is identical to the standard plane pair model. In spite of the significance test, the model is less robust to imperfectly filtered ground points or irregular point arrangements. In case of doubt the plane pair fit is used as a fall back strategy. Currently, the polynomial pair fit is only used if:

  • the curvature of the 2D approximations is low
  • the representative point of plane pair fit has a significant 2D offset to the 2D approximation
  • the horizontal σ of representative point is lower than the one from the plane pair fit
Fig. 2c: Polynomial patch pair fit

Individual point weights are used within patch pair fitting procedure. Points located in the middle of the patch hereby get a higher weight compared to points at the patch border. This applies both to the longitudinal and transversial direction. In both cases, a gaussian curve is used with weights of 1 in the middle and approximately 0.5 at the border. In addition to that, the height precision of the 3D points and the planimetric precision of the line approximations can be specified via &sigma a-priori. Basically, two values are expected (1st value=height precision of 3D points, 2nd value=planimetric precision of 2D line approximations). If only one value is provided, the planimetric precision is set to 3 times the height precision. Whereas the above mentioned parameters provide user control of the geometry pair fitting, the distance between to patch centers along the line can be defined via parameter overlap. Again, two values denoting the upper and a lower bound for the potential patch overlap are accepted. Internally, the algorithm adapts the actual overlap depending on the local surface shape. A large overlap is used in high curvature areas (e.g. at U-turns at the end of ditch, etc.) and a small overlap is employed if the line progression is rather straight.

As stated before, the course of the 3D structure line is modeled by intersecting local geometry primetives wherever possible. However, the planimetric uncertainty increases with decreasing intersection angles. Therefore it is possible to specify a critical intersection angle (minAngle). If the patch pair intersection angle is smaller than the specified value, the elevation is derived from a one-sided plane and the planimitric position of the 2D line approximation is preserved. If plane fitting succeeds for both the left and the right plane, the more horizontal plane is used for the heigth estimation. In case the output is requested in ODM file format, all result parameters of the geometry fit are stored for each patch center point. Please note, that no matter which output format is chosen for the final 3D structure lines, the ODM containing the above mentioned metadata is stored in any case. The most important parameters are:

  • standard deviation of geometry fit
  • used geometry model (i.e., plane pair, plane-cone, pair of polynomial cylinders, independent plane, one-sided plane
  • intersection angle
  • 3D-tangent direction of structure line at the patch reference point
  • 3D-surface normal direction vector of both the left and the right patch side

After modeling an entire line as a sequence of local patches of arbitrary length, a homogeneous point spacing along the 3D strucutre line is established by interpolating intermediate line vertices in approximately regular intervals (samplingDist). Cubic Bezier curves based on the 3D patch reference point positions and the corresponding 3D tangent directions of two adjacent patches are used for the line densification. Whereas generally all lines (cf. processId and ignoreId) above) are processed only those lines meeting a miminum user-defined length (minLength) are considered for output to suppresss overly short lines.

Line quality assessment

Once processing of a single line is completed, a quality assement is carried out. Quality indicators are provided as aggregated statistics for the entire line which might consist of severel (continuous) line parts (due to data voids in the 3D pint cloud, etc.). The quality level is described in school grade catagories (very good, good, moderate, sufficient, not sufficient). An additional class is provided for all lines which are topologically inconsistent (self intersections, etc.). The classification rule is as follows:

  1. Very good
    • Mean intersection angle > 1.5 times minAngle
    • Mean σ of geometry fit < a-priori height precision (sigmaApriori)
    • Maximum σ of geometry fit < 1.5 times a-priori height precision (sigmaApriori)
    • Number of line parts == 1 (i.e. continuous modeling of, no line breaks)
    • No one-sided or independent plane estimations
  2. Good
    • Mean intersection angle > 1.25 times minAngle
    • Mean σ of geometry fit < a-priori height precision (sigmaApriori)
    • Maximum σ of geometry fit < 2.0 times a-priori height precision (sigmaApriori)
    • Number of line parts < line length / 5 times minLength (i.e. few line parts)
    • Less than 5% one-sided or independent plane estimations
  3. Moderate
    • Mean intersection angle > 1.0 times minAngle
    • Mean σ of geometry fit < 1.5 times a-priori height precision (sigmaApriori)
    • Maximum σ of geometry fit < 2.5 times a-priori height precision (sigmaApriori)
    • Number of line parts < line length / 5 times minLength (i.e. few line parts)
    • Less than 10% one-sided or independent plane estimations
  4. Sufficient
    • Mean intersection angle > 0.75 times minAngle
    • Mean σ of geometry fit < 2.0 times a-priori height precision (sigmaApriori)
    • Maximum σ of geometry fit < 2.5 times a-priori height precision (sigmaApriori)
    • Number of line parts < line length / 3 times minLength (i.e. many line parts)
    • Less than 20% one-sided or independent plane estimations
  5. Unsufficient

    • otherwise

    Please note that, as mentioned before, more detailed metadata (on patch level) are available in the ODM output file.

Final output

The structure lines are provided in (i) ODM, and (iia) ESRI shapefile, or (iib)(Binary) Winput file format. In case of (Binary) Winput a separate output file is created for each catagory (very good, good, etc.). Furthermore, Winput does not provide the capability to store additional attributes. For ESRI shapefiles, a single file is created and the line id as well as the quality level (_origId, _quality) are provided as attributes in the corresponding database (dbf) file. Furthermore, a secondary output file containing the individual line segments is provided. For each segment the following attributes are available: (i) _id, (ii) _segmentId, (iii) _origId, and (iv) _curvature. Whereas _id constitutes a countinuous counter, _segmentId denotes the index of a segement within a certain line, _origId the id of the corresponding 3D structure line and, finally, _curvature a concave/convex indicator.

Parameter description

-inFileODM input file
Type: opals::Vector<opals::Path>
Remarks: mandatory
The point cloud defined by the specified ODM input files represent the central data pool for extracting and modelling structure lines.
-approxFileapproximative structure lines
Type: opals::Path
Remarks: mandatory
Based on the given 2D lines in approximation file the extraction and modelling of the structure lines is performed.
-outFileoutput file of modelled structure lines
Type: opals::Path
Remarks: mandatory
Based on the given lines in approximation file the extraction and modelling of the structure lines is performed. Estimation rule: The current path and the name (body) of approxFile postfixed by '_modeled'
-oFormatoutput format [odm, wnp, bwnp, shape]
Type: opals::String
Remarks: default=shape
Beside an ODM, which is always created, the modelled lines can be additionally exported as shape or winput file
-filtersub-select input ODM
Type: opals::String
Remarks: optional
A filter string in EBNF syntax as described in section 'Filter syntax' can be passed to restrict the inputpoint cloud ODM (e.g. to consider last echoes only)
-processIdprocess specific line ids only
Type: opals::Vector<opals::String>
Remarks: optional
The value can either be (multiple) line ids or a single textfile containing multiple line ids which will be processed. All other lines are skipped.
-ignoreIdignore specific line ids
Type: opals::Vector<opals::String>
Remarks: optional
The value can either be (multiple) line ids or a single textfile containing multiple line ids. Lines with those ids will be skipped during processing.
-patchLengthmin and max patch length
Type: opals::Vector<double>
Remarks: default=5 15
Sets the minimum and maximum patch length for modelling a single patch. Based on the curvature of the approximateline the patch length is varied between the given minimum (straight line) and maximum (high curvature) value. If only onevalue is provided, the curvature based variation is omitted.
-patchWidthleft and right patch width
Type: opals::Vector<double>
Remarks: default=5
Sets left and right patch width for modelling a single patch. If only one value is given, it is used for left andand right patch width. Hence, the width of the actual patch rectangle is left+right width.
-overlapoverlap of adjacent patches
Type: opals::Vector<double>
Remarks: default=0.15 0.75
Sets the minimum and maximum overlap between adjacent patches. Based on the curvature of the approximateline the patch overlap is varied between the given minimum (straight line) and maximum (high curvature) value. If only onevalue is provided, the curvature based variation is omitted.
-minAngleminimum intersection angle [deg]
Type: double
Remarks: default=7
The minimum angle between left and right surface is mainly used for quality assessment of the final line.
-minLengthminimum length of final line
Type: double
Remarks: default=0
If modeled line is shorter than the given threshold, it will be ignored in the final export.
-samplingDistsampling interval of final line
Type: double
Remarks: default=1
The vertices of modeled line will be evenly spaced with the given sampling interval.
-pointCountmin und max point count for patches
Type: opals::Vector<uint32>
Remarks: default=10 0
The first value is the minimum number of points required for doing a model fit left or right of the approximation line. If a second value is specified (or a value other than 0), it is used as maximum point count per patch resulting in local patch length reductions.
-sigmaAprioripoint height and line approximation sigma
Type: opals::Vector<double>
Remarks: default=0.15
The first value defines the sigma of the height observations for all model fits. The second value sets the 2D sigma of line approximations. If no second value is specified the height sigma is multiplied by 3 and used as 2D sigma.
-spikeRemovalspike removal flag
Type: bool
Remarks: default=1
If enabled the process removes 2d spikes in the final modelled line. In case of awkward point distribution or weakly defined structure lines spikes may occur if different models are used for neighboring patches.
-debugOutFileprefix for patch modelling results
Type: opals::Path
Remarks: optional
If specified debug winput files are created for visually analysing the patch results. Multiple output files will be created so only a file prefix needs to be specified (The winput code translates to the used adjustment model: 50 -> plane pair, 40 -> plane cone, 12 -> polynomial pair, 60 -> single plane)

Examples

The data used in the following example can be found in the $OPALS_ROOT/demo/ directory. As a prerequisite the data of (flyover.laz) are imported into an OPALS data manager, and 2D line approximations are derived following the procedure detailed in the examples section of Module LineTopology.

opalsImport -inFile flyover.laz -tileSize 100 -coord_ref_sys EPSG:5650
opalsGrid -inFile flyover.odm -outFile flyover_dtm.tif -interpolation movingPlanes -neighbours 8 -searchRadius 2 -selMode quadrant
opalsGridFeature -feature slpDeg -inFile flyover_dtm.tif
opalsEdgeDetect -inFile flyover_dtm_slpDeg.tif -outFile flyover_dtm_canny.tif -sigmaSmooth 1.25 -threshold 1 4
opalsVectorize -inFile flyover_dtm_canny.tif -outFile flyover_dtm_canny.shp
opalsLineTopology -inf flyover_dtm_canny.shp -outf flyover_dtm_canny_step1.shp -cfg linetopology_step1.cfg
opalsLineTopology -inf flyover_dtm_canny_step1.shp -outf flyover_dtm_canny_step2.shp -cfg linetopology_step2.cfg
opalsLineTopology -inf flyover_dtm_canny_step2.shp -outf flyover_dtm_canny_step3.shp -cfg linetopology_step3.cfg

The first example reads the paprameters from the configuration file (lineModeler.cfg) provided in the demo directory.

oFormat=shape
filter=class[ground]
patchLength=4.0 15.0
patchWidth=2.5 2.5
overlap=0.15 0.75
pointCount=10 350
minAngle=7.0
minLength=10.0
samplingDist=1.0
sigmaApriori=0.10 0.25

The 3d structure lines can then by derived simply by passing the path and name of the configuration file as well as the path names of the point cloud (inFile), of the 2d line approximations (approxFile) and of the resulting 3D structure lines (outFile).

opalsLineModeler -cfgFile lineModeler.cfg -inFile flyover.odm -approxFile flyover_dtm_canny_step3.shp -outFile flyover_3d_breaklines.shp

The resulting 3D structure lines are plotted in Figure 1 together with the original 3D point cloud.

Fig. 3: 3D structure lines (blue) derived from 2D line approximations and 3D point cloud

A quality measure for the resulting lines is the distance to the DTM surface. It can be calculated for each node of the structure lines. First, the nodes have to be exported as (xyz)-Points, and then imported to an odm again:

opalsExport -inFile flyover_3d_breaklines.odm -outFile flyover_breakl_points.xyz
opalsImport -inFile flyover_breakl_points.xyz -outFile flyover_breakl_points.odm

Using opalsAddInfo, the distances to the DTM can be calculated:

opalsAddInfo -inFile flyover_points.odm -gridFile flyover_dtm.tif -attribute "_dist2DTM(float)=r[0]-Z" -resampling bilinear

With a custom format definition file (addons\formatdef\lineModelerShape.xml), the attribute can be exported to a shapefile for further investigation:

opalsExport -inf flyover_points.odm -oformat $OPALS_ROOT\addons\formatdef\lineModelerShape.xml -outf flyover_dist2dtm.shp

References

G. Mandlburger, J. Otepka, C. Briese, W. Mücke, G. Summer, N. Pfeifer, S. Baltrusch, C. Dorn, H. Brockmann: Automatische Ableitung von Strukturlinien aus 3D-Punktwolken. in: "Dreiländertagung der SGPF, DGPF und OVG: Lösungen für eine Welt im Wandel", T. Kersten (ed.); Publikationen der Deutschen Gesellschaft für Photogrammetrie, Fernerkundung und Geoinformation e.V., Band 25 (2016), ISSN: 0942-2870; 131 - 142.

Otepka, J., Mandlburger, G., Briese, C. und Pfeifer, N.: Fortschritte bei der automatischen Ableitung von Strukturlinien aus 3D-Punktwolken. 19. Internationale Geodätische Woche Obergurgl 2017, Hanke, K. & Weinold, Th. (Hrsg.) Herbert Wichmann Verlag, VDE VERLAG GMBH, Berlin/Offenbach. ISBN 978-3-87907-xxx

Author
jo, gm
Date
1.3.2017
@ quadrant
quadrant-wise nn selection, ie. nn per quadrant, then 2nd nn per quadrant, ...
opalsVectorize is the executable file of Module Vectorize
Definition: ModuleExecutables.hpp:243
@ approxFile
File containing approximation info.
@ gridFile
grid model file
@ movingPlanes
moving (tilted) plane interpolation
@ coord_ref_sys
default coordinate reference system (EPSG Code, WKT string or PRJ-File)
opalsAddInfo is the executable file of Module AddInfo
Definition: ModuleExecutables.hpp:8
@ bilinear
bi-linear interpolation
@ cfgFile
configuration file
opalsImport is the executable file of Module Import
Definition: ModuleExecutables.hpp:113
opalsEdgeDetect is the executable file of Module EdgeDetect
Definition: ModuleExecutables.hpp:68
@ threshold
threshold (opalsEdgeDetect)
opalsExport is the executable file of Module Export
Definition: ModuleExecutables.hpp:73
opalsLineModeler is the executable file of Module LineModeler
Definition: ModuleExecutables.hpp:133
opalsGrid is the executable file of Module Grid
Definition: ModuleExecutables.hpp:93
opalsLineTopology is the executable file of Module LineTopology
Definition: ModuleExecutables.hpp:138
@ slpDeg
steepest slope [deg]
Definition: GridFeature.hpp:18
@ sigmaSmooth
gaussian smoothing sigma (opalsEdgeDetect)
@ feature
Use a statistic feature of the boundary gap points for filling.
opalsGridFeature is the executable file of Module GridFeature
Definition: ModuleExecutables.hpp:98