Aim of module

The aim of the module TerrainFilter is to classify a point cloud (ODM) into terrain and off-terrain points based on the hierarchical robust interpolation (Kraus and Pfeifer, 1998; Pfeifer et al., 2001; Pfeifer and Mandlburger, 2018) approach using robust surface interpolation with kriging. The classification result is stored in the point attribute Classification. To visualize the results, use Module View, or create raster DTM with Module Grid and a shaded relief map with Module Shade (cf. terrainfilter_examples examples).

General description

Module TerrainFilter implements the hierachic robust interpolation method described in Pfeifer et al. (2001) and Pfeifer and Mandlburger (2018). The core of the approach is an iterative robust interpolation procedure (Kraus and Pfeifer, 1998) based on linear prediction (kriging), which can be used to separate terrain from off-terrain points. The method starts with an initial approximation of the surface using all points (vegetation, powerlines, buildings, etc., and ground). Next, the residuals, i.e. the oriented distances from this approximate surface to the measured points, are computed. Each point is given a weight according to its distance from the surface, priorizing points below over points above the surface. The surface is then recomputed under the consideration of the weights. A point with a high weight will attract the surface, resulting in a small residual at this point, whereas a point that has been assigned a low weight will have little influence on the run of the (subsequent) surface. This process of weight iteration is perpetuated until all gross errors are eliminated (a stable situation) or a maximum number of iterations is reached. Finally, a point can be classified as off-terrain, if its oriented distance is above or below certain values. Figure 1 illustrates the iterative robust interpolation approach.

Figure 1: Iterative robust interpolation; (a) initial surface + weight function; (b) surface after 1st iteration; (c) final surface and classified terrain and off-terrain points

In order to classify large unpenetrable objects like houses and factory buildings, robust interpolation is embedded in hierarchic approach (Pfeifer et al., 2001; Pfeifer and Mandlburger, 2018). The (spatial) hierarchy is implemented in Module TerrainFilter by specifying the smallest unit via parameter robustInterpolation.gridSize and the number of hierarchy levels via parameter robustInterpolation.pyramidLevels. The edge lengths of a (regular) grid cell in the subsequent hierachy level is twice as large as the grid cell of the current level. Thus, the grid size of the coarsest hierachy level $gs_{n-1}$ computes to $2^{n-1}$, with $n$=robustInterpolation.pyramidLevels. Thus, choose the number of pyramid/hierachy levels such that the largest building is smaller than $gs_{n-1}$. For robustInterpolation.pyramidLevels=4, Figure 2 illustrates the hierarchic grid structure.

Figure 2: Expemplary structure of hierarchic grid for 4 pyramid levels

Hierarchic filtering starts with the coarsest pyramid level, and in each but the finest level the following steps are performed (cf. Figure 3):

First, a representative point is determined in each grid cell from all active points using a certain feature (min, n:min, certain quantile; parameter robustInterpolation.feature). In the coarsest level, all points passing the user defined filter (parameter filter) are used. In all subsequent hierarchy levels, points outside a certain height tolerance band might be descarded (cf. below).
Robust interpolation is performed with the selected points, resulting in an approximative surface for each hierarchy level.
A height tolerance band around this surface is constructed. The upper thresholds (i.e. distance above the surface) can be defined individually for each hierachy level via parameter robustInterpolation.filterThresholds. In general, larger values are used for the coarser levels (i.e. mean building height in the last level) and gradually smaller thresholds should be used for the finer levels. The height band does not have to symmetric around the surface, but the lower bound is typically larger than the upper bound. The asymmetry can by specified specified via robustInterpolation.lowerThresholdScale (default: -1.5 * upper threshold).
Finally, all points outside the height tolerance are deactivated for the next hierachy level (but can become active again in the subsequent levels).
The entire procedure is repeated for the next finer pyramid level with all points within the height tolerance band.

Figure 3: Sequence of processing steps of a single hierarchy/pyramid level; (a) Creation of data pyramid and selection of representative points per grid cell; (b) Iterative robust interpolation with all representative points; (c) Construction of height tolerance band; (d) Selection of active points for next hierarchy level

Please note that in the last (=finest) pyramid level, the tolearance band (parameters robustInterpolation.filterThresholds, robustInterpolation.lowerThresholdScale) are used for the final classification of the points. The classification attribute of the points within the tolerance band is changed to 2 (i.e., ground) while points the classification of all points outside the band remains unchanged.

For a detailed explanation of all remaining parameters, please refer to Section Parameter description.

Parameter description

-inFileinput ODM file

Type: opals::Path
Remarks: mandatory
The path to the opals datamanager whose point data are being classified

-methodterrain filter method

Type: opals::TerrainFilterMethod
Remarks: default=robustInterpolation

Possible values:  
  robustInterpolation ... terrain points filter via (hierarchical) robust surface interpolation

Currently only hierarchical robust surface interpolation is supported.

-limit2D clipping window

Type: opals::GridLimit
Remarks: optional
If no user defined limits are specified or -limit is even skipped, the OPALS Data Manager (ODM) extents are used. Rounding is enabled by default in case the limits are derived from the ODM and can be enforced for user defined limits (keyword: round).
GridLimit: Defines xy-limits for output grid/raster datasets. The limit coordinates specified by left/lower/right/upper either refer to pixel centers (default) or corners, and may optionally get rounded to multiples of the grid size.
Syntax: ['center'|'corner']['round']['(' left lower right upper ')']

-filter(tree of )filter(s) to select data points/lines

Type: opals::String
Remarks: optional
If a filter string is specified, only those points and lines passing the filter take place in the actual robust interpolation. Hovever, all other points (within the specified limit window) are also classified into terrain and off-terrain (i.e. last step of processing pipeline).

-tempDirectorypath for intermediate files

Type: opals::Path
Remarks: estimable
Directory path for storing intermediate files).
Estimation rule: Current working directory + '/tempTerrainFilter'

-deleteTempDatadelete intermediate data files

Type: bool
Remarks: default=0
Temporary data are stored within the directory -tempDirectory. If true, the temporary directory is delete at the end of the run.

-debugOutFiledebug terrain points file

Type: opals::Path
Remarks: optional
Outputs a xyz file containing all points classified as ground

-writeFilterInfowrite filter information to the odm

Type: bool
Remarks: default=0
Writes additional information about the filter process to user-defined attributes of the odm.

-sigmaAprioriIGroup: A-priori sigma options

.bulkPointsheight accuracy of bulk points

Type: opals::String
Remarks: default=0.10
Either a constant accuracy representing all bulk points or a user-defined formula using the generic filter syntax may be applied. It is important thatthe a specified formula provides realistic accuracy values, otherwise therobust interpolation concept will not work properly.

.keyPointsheight accuracy of key points

Type: opals::String
Remarks: default=0.05
Detailed description cf. above (sigmaApriori.bulkPoints).

.formLinesheight accuracy of form line vertices

Type: opals::String
Remarks: default=0.03
Detailed description cf. above (sigmaApriori.bulkPoints).

.breakLinesheight accuracy of break line vertices

Type: opals::String
Remarks: default=0.03
Detailed description cf. above (sigmaApriori.bulkPoints).

-robustInterpolationIGroup: Robust interpolation options

.gridSizebasic grid cell size

Type: double
Remarks: optional
Hierarchical robust interpolation works from coarse to fine in multiple resolutions (pyramid levels). If set, the specofied value denotes the grid width of the finest resolution (i.e. 1st pyramid level). If not set, the respective grid width is estimated based on the point density of the data.

.pyramidLevelsnumber of data pyramid levels

Type: uint32
Remarks: default=4
The hierarchical robust interpolation strategy works from coarse to fine. Whereas the resolution of the finest surface grid model is derived automatically by analyzing the point density of the dataset, the coarser levels are derived by repeatedly applying a 2x2 kernel.

.featurefeature for point selection in pyramid levels

Type: opals::Vector<DM::StatFeature>
Remarks: optional
For each pyramid level a different statistical feature used for selecting a representative point within a grid cell can be specified. In general, the lowest points within grid cells(min) represent good starting points for deriving the terrain surface with robust interpolation. If long range outliers are present in the dataset, a more robust feature (e.g. quantile:0.05 (=default) or nmin) may be advantageous. Specify 'null' fo suppress thinning for a specific level (recommended for first (i.e. full resolution ) level only.

.filterThresholdsupper filter threshold in pyramid levels

Type: opals::Vector<double>
Remarks: default=0.2 0.5 1 3
Either specify the upper filter treshold for each pyramid level or specify the upper filter threshold for the first (max) and the last pyramid level (min). All other upper thresholds will then be computed by linear interpolation. The lower filter threshold is computed for a level i, based on the filterThreshold with the parameter 'lowerThresholdScale' as filterThreshold(i)*lowerThresholdScale.

.lowerThresholdScalescale for the lower threshold

Type: double
Remarks: default=-1.5
The lower filter threshold is computed for a level i, based on the parameter 'filterThreshold' as filterThreshold(i)*lowerThresholdScale.

.maxSigmamaximum allowed sigma of interpolation

Type: double
Remarks: default=0.5
If the standard deviation of a unit observation of the robust interpolation (i.e. sigma_0 ) exceeds the specified value, the classification of the affected points remains unchanged.

.penetrationestimated penetration rate [%]

Type: uint32
Remarks: default=20
The laser signal is partly or totally reflected at semi-transparent objects like vegetation. Therefore, not all laser pulses reach the ground. The penetration rate is used in the course of the robust interpolation for apriori estimating a reasonable initial course of the local surfaces

.maxItermaximum number of iterations

Type: uint32
Remarks: default=10
Robust interpolation is performed iteratively. In each iteration the individual point weights are adapted based on the residuals (i.e. vertical distance between point and intermediate surface. The process of surface interpolation and re-weighting is repeated until the surface changes are below a threshold or the maximum number of iterations is reached.

.tileSizetile (block) size

Type: uint32
Remarks: default=128
The internal processing is organized in rectangular tiles. The size of these processing unitscan be specified allowing control of the memory consumption of the module. Please note, that the tileSize value is interpreted in units of the first pyramid level grid width rather than in metric units. The value should be larger than the largest building of the dataset.

.robustWFAdpationadaption of the robust weight function

Type: opals::String
Remarks: default=adapting
The module supports different adaption functions

-classifyOverlapclassify points in overlap multiple times

Type: bool
Remarks: default=0
Each CU has an overlap with its neighbors (netto vs. brutto size). If this flag is enabled, the all points in the brutto CU are classified, meaning that points in the overlapping area get potentially classified multiple times.

Examples

The data used in the following examples can be found in the $OPALS_ROOT/demo/ directory. The point cloud of fullwave2.fwf is imported into an OPALS data manager using the following command:

opalsImport -inFile fullwave2.fwf

Example 1

The first example classifies the 3D ALS point cloud of a forest scene. Since no buildings are available, it is sufficient to use only 3 pyramid levels with grid sizes of 1m, 2m and 4m, respectively. In the final level all points within a band of -0.75 ot +0.5m are accepted as ground points.

opalsTerrainFilter -inf fullwave2.odm -filter echo[last] -pyramidLevels 3 -filterThresholds 0.5 1 2 -gridSize 1

The classified point cloud is displayed in Figure 4.

Figure 4: Z-colored fullwave2.odm dataset and the ground points (Attribute: Classification)

To compute a raster DTM and the shading the following commands can be used:

opalsGrid -inf fullwave2.odm -outFile fullwave2.tif -interpolation movingPlanes -searchRad 2 -filter "class[ground]"

opalsShade -inf fullwave2.tif

Example 2

The second example uses the demo dataset strip11.laz, i.e. a residential area with 2 storey houses densly packed along the main road. In order to filter the buildings, 5 pyramid levels are necessary, resulting in a grid size of the first/coarsest level of $2^4=16m$ (using the default grid size of 1m). In the coarsest level an upper tolerance of 4m exludes all roof points for the subsequent levels. In the final level, a small threshold of 0.2m provides prevents low vegetation (grass, bushes) from being classified as ground. Filtering can be achieved with the following command sequence:

The point cloud of strip11.laz is imported into an OPALS data manager using the following command:

opalsImport -inFile strip11.laz

To reset an existing classification on the dataset use the following command:

opalsAddInfo -inf strip11.odm -attr "Classification=invalid"

opalsTerrainFilter -inf strip11.odm -filter echo[last] -pyramidLevels 5 -filterThresh 0.2 0.5 1 2.5 4

Figure 5: Z-colored strip11.odm dataset and the ground points (Attribute: Classification)

References

Kraus, K. and Pfeifer, N., 1998. Determination of terrain models in wooded areas with airborne laser scanner data. ISPRS Journal of Photogrammetry and Remote Sensing 53, pp. 193-203.
Pfeifer, N., Stadler, P. and Briese, C., 2001. Derivation of digital terrain models in the SCOP++ environment. In: Proceedings of OEEPE Workshop on Airborne Laserscanning and Interferometric SAR for Detailed Digital Terrain Models, Stockholm, Sweden.
N. Pfeifer, G. Mandlburger: Lidar data filtering and Digital Terrain Model generation. in: Topographic Laser Ranging and Scanning - Principles and Processing, Second Edition, J. Shan, C. Toth (ed.); CRC Press, Boca Raton, 2018, ISBN: 9781498772273, 349 - 378.

Author: gm, mpoechtr

Date: 04.06.2022

Table of Contents