oht 
The 
| to 
be 
one 
lere 
| to 
) tO 
Iter 
set- 
een 
und 
‘the 
e of 
hich 
nof 
'eral 
f all 
(1) 
veen 
Alter 
(4) 
George Vosselman 
€ =min|h, AD, s (ap, p,))] (5) 
PEA 
Now the relationship between the erosion and the defined ground point filter can be expressed. If h, = e, , then 
Vp,€ A:h,, Sh, * Ah, (d(p,.p)) (6) 
Thus, the set of ground points can be defined by 
DEM -|p, € AIh, xe, ] (7) 
. In words: a point is classified as a ground point if its height does not exceed the height of the eroded surface. Filtering 
laser altimetry data can therefore be done with two elementary operations: a morphological erosion with the kernel 
function as defined in (4) and a test on the difference between the original height and the eroded height of a point. 
2.3 Implementation issues 
Morphological operators are available in most image processing packages. These packages usually require the data to 
be organised in a grid. Interpolating the raw laser altimetry data to a grid, however, causes a significant loss of 
information. In particular when heights are interpolated between ground points and points on vegetation or buildings, 
the height differences in the interpolated data will be reduced. Therefore, it becomes more difficult to make a correct 
classification. Although it is computationally more expensive, we prefer to work with the original, irregularly 
distributed point data. 
According to the filter definition in equation (1), the height of a point needs to be compared with the heights of all other 
points. In most cases this is not necessary. For example, if it is known that the height difference in the terrain is no more 
than 10 m and that Ah, (100 m) = 10 m, only points within 100 m of the current point need to be considered. In order 
to determine the points within some distance, we organise the points in a Delaunay triangulation. 
Instead of verifying that a point fulfils the condition in equation (6), one can also check whether the height of a point 
results into the rejection of a point in its neighbourhood. This approach has been taken in the implementation. First the 
height of the direct neighbours in the triangulation are verified. Similar to a region growing algorithm, in the next step 
the neighbours of these neighbours are verified. The neighbourhood is grown until there are no more adjacent points 
within the maximum distance that needs to be considered. To speed up the filtering process, the neighbours of a node 
are only verified if the node itself is rejected. This implies that points in the DEM set may not strictly comply with the 
definition in equation (7). The number of errors introduced by this heuristic is very small. The processing time, 
however, reduces with a very large factor (10-100). 
3 DERIVATION OF FILTER KERNELS 
The filter function should incorporate the knowledge about the height differences in the terrain. In this section we 
describe three ways to derive a filter function. The first one uses simple generic knowledge about the terrain shape and 
the precision of the height measurements. The other two make use of a training data set. 
3.1 Synthetic function 
Suppose that, for some area, one would know that the slopes in the terrain are not steeper than, say, 30%. If the 
measurements are free of errors, the filter function could be defined as: 
Ah (d)=03d (8) 
In the usual case of noisy measurements, one would like to add a confidence interval. If one would allow that 5% of the 
terrain points with a standard deviation 6 may be rejected, the filter function becomes 
Ah ux (0) 2 0.3d 11.654260 (9) 
  
International Archives of Photogrammetry and Remote Sensing. Vol, XXXIII, Part B3. Amsterdam 2000. 937