'LE PASS
jects belonging
sensors makes
orithm has been
Digital Terrain
ot the ground is
gorithm moves
ated ground and
| ground surface
al DTM will be
roceed relevant
s with different
ze over a DEM.
eir belonging to
it in the applied
round, whereas
ind points.
morphological
elt, 1995). The
h a mask in as-
orresponds to a
rived using this
non-ground ob-
ach may be ex-
y means of air-
, Z), the dilation
ined as:
in (Guy: 201)
p)Ew
(z,y,2) and
described dual-
t al., 2003) ina
rs increase (ex-
ies of recursive
und if the ele-
ation k and the
that depends on
1998) is based
erage surface to
given a weight
eter of a weight
International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol XXXV, Part B3. Istanbul 2004
function. The surface is then recomputed under the consideration
of the weights. Intuitively, it is assumed that terrain points are
more likely to have negative residuals, whereas vegetation (build-
ings) points are more likely to have positive residuals. During
these iterations a classification is performed. If an oriented dis-
tance is above a certain value, the point is classified as off-terrain
point and eliminated completely from this surface interpolation.
Lee (Lee and Younan, 2003) modified the previous method by im-
plementing an adaptive prediction technique for extracting DTMs
of the ground surface underlying vegetation. According to the au-
thors, this technique offers, in general, a better tracking capability
in the extraction of bare Earth models with steep slopes and large
variability.
Surface based
An other approach was introduced in (Axelsson, 1999) based on
the connection of a surface from below the point cloud. This sur-
face is connected to the ground points using different criteria such
as the Minimum Length Description (MDL), constrained spline
functions or snakes. All criteria are meant to manage the possible
shapes and hence the fluctuations of the resulting surface. The
active shape models were first applied to Lidar data in (Elmqvist
et al., 2001) and (Elmqvist, 2002). Raw data are first re-sampled
over a regular grid, the ground surface is then estimated with the
minimization of a defined energy which depends on the inter-
nal behavior of the surface, on a data term and on other external
forces. According to the author, this algorithm is very robust, and
it works on data of different types of terrain.
An other method which is continuously adaptive to terrain sur-
face variations has been developed (Sohn and Dowman, 2002). It
aims to recursively divide the LIDAR data into a set of piecewise
planar surface models in order to underly terrain slope variations
regularized into homogeneous plane terrain. The authors used
a downward divide-and-conquer triangulation to run in the point
cloud.
Geometry
The slope-based filter uses the slope of the line between any two
points in a point set as the criteria for classifying ground points
(Vosselman, 2000). If the slope exceeds a certain threshold then
the highest point is assumed to belong to an object. This filter
was modified so that the threshold should vary with respect to the
slope terrain (Sithole, 2001).
3 DESCRIPTION OF THE ALGORITHM
The algorithm is based on a bipartite voting process. A laser
point will be labeled several times either as ground or non-ground
points until the most represented label be affected to the final clas-
sified point. Following the propagated direction (section 3.1), an
estimation of the ground is performed (section 3.2). This prime
DTM is then refined by using an energy minimization algorithm
so that the final DTM should be as accurate as we may expect
from laser data.
3.1 Propagation
The propagation mechanism consists of moving onto a regular
geocoded grid. Starting from the cell whereupon the lowest laser
point is included, the algorithm explores his 4 non-visited grid
neighbors (4-connexity). It then extracts the corresponding laser
neighbors V (equation 2) and insert both their average altitude
and their position in a sorted (ascending order) container struc-
ture.
Ik
/ = = [conter = Tul < C 9
j Ux Uk 1 [o zt ra < e } Q)
2k kEN
where ly, is a laser point, (X center; Ycenter ) is the planimetric cen-
ter of the neighborhood and C is a constant. The algorithm prop-
agates itself toward the lower feature of this structure. Figure 1
sketches the behavior of the propagation: black cells have already
been visited whereas gray ones are potential candidates.
Figure 1: Aspect of the propagation on a geocoded grid. Black
cells have already been visited whereas gray ones are potential
candidates
3.2 Segregating bare/non-bare earth points
How a laser point is temporary classified as a terrain point? An
altimetric difference 1s calculated between the z component of
the laser point and the estimated terrain elevation at this place
(round local). This estimation is a mean of the 2096 lowest
laser points of the neighborhood. If the difference is less than a
fixed threshold (say 50cm), the point is considered to belong to
the terrain. A new estimation of the terrain height is then com-
puted (mean of points classified as ground) taking into account
the new ground laser point. If the difference is larger, the point
is classified as a non-ground point. This calculation is performed
until all laser points belonging to Vi,; ((?, j) are the coordinates
of the geocoded grid) be processed. A prime DTM (denoted
Sin in this paper) is then filled at (7, 7) with the ground value
Zorou nd local:
The algorithm has its own error self detector which will un-
derline both erroneous ground estimations and 3D-point mis-
classification. This detector is based on the comparison between
the average local elevation of the ground 2 round local OF Sui
(3 x 3 window size without the central cell, see Figure 2(a)) and
the above calculated Zground tocat. A point classified as ground
will be detected if the local slope is larger than arctan ak, where
R is the resolution of S;n and Ah is the altitude difference. We
therefore apply a linear correction (equation 3) to the current
ground elevation in order to take into account the real local slope.
^
Zgroundiocal € aZground local + (1 QYZ around local (3)
where a a constant chosen depending on the respective weight
we want to grant either t0 Z ground local OF 10 Zoround local:
Considering the overlapping structure of our neighborhoods,
laser points are classified several times (exactly (£)7). As we
can see on Figure 3, V;,; and V;i41,; have a large number of com-
mon points (empty circles). As a result, for each neighborhood
extraction, points will be labeled following local criteria. At the
end of the propagation, a laser point will have been labeled p
times as ground and m times as non-ground. We then affect the
final label corresponding to max(n, p), which is the most repre-
sentative vote.
