spground

Classifies points as ground and non-ground, and can either (-groundonly) output ground points only, (-nongroundonly) output non-ground points only, or (default) output both, with non-ground points translated upward so they can be easily distinguished by downstream tools. (This is a hack, but it is less expensive in terms of storage than keeping a type with each point, and thus we have used it extensively in these prototype tools.)

Example Usage:

Usage:

   spground.exe -i input.spb -o split_points.spb -tmax 30

Other parameters:

The algorithm is described below. It depends on three parameters that can be adjusted by command line flags

 -slopelimit 0.5 
 -tmax 30
 -tmin 2 

Details:

We have used a filtering algorithm to classify the ground points from those belonging to non-ground objects (buildings, cars, traffic lights, power lines, other man-made objects). The algorithm was initially developed by Lindenberger at al [1] and used by the company TopScan [2]. The most detailed description we could find is from the paper of Petzold at al. [3]. Then we have modified and adapted the algorithm, so that it can handle large LIDAR datasets without keeping all the data in the memory.

Description of the original algorithm:

The algorithm works by iteratively finding and discarding a set of points which are situated above the ground.

• First, the entire data set is scanned by a moving window and the points with lowest z-value in the window are captured. The collection of all these points represents a rough estimation of the DEM.
• Then we can classify all points above certain threshold from the DEM to be non-ground. All these points are removed from the dataset.

We repeat the same process with smaller window size (and smaller threshold) and this gives us a more detailed representation of the ground. Finally, when we reach a small enough size for the window, we can assume that we have discarded all non-ground points, keep the rest and mark them as ground points.

Adapting the algorithm to streaming datasets:

The algorithm, as defined above, imposes some conditions:

• to have quick access to a given subregion (the window). This usually means that we have all the data read and stored somewhere.
• to store all the lowest points. We need them to filter out the non-ground points at each step. The number of these points is proportional to the total size of the dataset, and for small window and flatty terrain, this number could get quite large.

Both requirements are in conflict with the streaming approach, where we want to keep the number points in memory to a minimum (and this number is practically independent of the dataset size).

First version:

We divided the input grid in NxN tiles (each tile around 100x100 meters). We read points from the stream along with finalization tags, and whenever we have read all the points for a given tile, we run the original algorithm on it. The output of that is a subset of the input tile data - these points which are detected as ground. We just output them preserving the finalization tags between them, which would preserve the finalized data format for whoever is is using this output. One simplification we make to the original algorithm is that the window size does not move continuously over the data set, but instead a number of fixed locations is taken (either a grid with the size of the window). We tried with a finer grid - 1/2 the size of the window, but the redult doesn't differ much. This works pretty well both in terms of efficiency and correctness of the classification. The obvious problem is that these tiles are making the output discontinuous at boundaries. Moreover, we want to keep the tile size small, so that memory consumption is kept at minimum - it is about 3 times the size of the initial window size. This leads to a relatively big portion of points which are on the boundary and are being evaluated without having sufficient neighborhood information.

Second version:

Our goal here is to remove the boundary artifacts from the first version and make the solution entirely 'streamed'. To do this, we make a imaginary division of the dataset in a grid with cell size w0 (the initial size of the window). We'll also take divisions with cell size w1, w2, ... (subsequent window sizes). For our implementation the window sizes decrease two times at each iteration, and this guarantees that the grids fit nicely into each other. What do we do with these grids:

• first, we want to find the lowest point at each cell in the grid with size w0.
• the difference with the original algorithm is that we don't process all cells, before going to the lower level (grid with size w1). Instead we can process a cell at w0 as soon as we have the lowest points of its 8 neighboring cells: then we just compute which points to discard and assign the preserved points into their respective w1 cells. Again, these w1 cells can be processed as soon as we have the lowest points at their 8 neighbors at w1 level.
• thus at any given time, each point active in the memory would belong to a level w0,... wk. Once we assign points to the lowest level, we can output them as ground points with the corresponding finalization tag.

This approach effectively processes the dataset looking at it as a whole, without the need to divide it into patches and glue them later. At the same time it keeps in memory only the active cells at each level. Thus the memory consumption is proportional to the one on the first version.

Results:

We ran the implementations on the 5x5 square:

First version:

• window size: w0 = 20.83m, w4 = 1.30m.
• memory used: 30mb.
• run time: 77sec.

Second versions:

• window size: w0 = 31.25m, w4 = 1.95m.
• memory used: 57mb.
• run time: 160sec.

• window size: w0 = 15.63m, w3 = 1.95m.
• memory used: 36mb.
• run time: 116sec.

Limitations:

One limitation of the original algorithm is that discontinuities on the ground are not detected well. The reason is that we use (linear) interpolation between the lowest points, and the DEM estimation would be some more or less smooth surface. In the data set, there was a place where there were two or three ground levels, possibly some construction site. Around the unevenness the points were classified as non-ground. To fix this, we would have to address the problem formulation - i.e. how exactly do we want to define ground? Then there were problems of small spotty areas being classified incorrectly. This could probably be fixed by tuning the parameters to the data set (initial window size, the thresholds, number of iterations). Since we might not be able to completely remove this 'noise' in the output, a more robust solution would be to post-process it.

1. J. Lindenberger, 1993. Laser-Profilmessungen zur topographischen Gelandeaufnahme. PhD dissertation, Institute for Photogrammetry, Stuttgart University.
2. http://www.topscan.de.
3. B. Petzold, P. Reiss, W. Stossel - Laser scanning - surveying and mapping agencies are using a new technique for the derivation of digital terrain models. ISPRS Journal of Photogrammetry & Remote Sensing 54(1999)95–104

-- JackSnoeyink - 12 Dec 2008
to top