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