B3. Istanbul 2004
identified and se-
ılgorithm on state-
n section 4.
rresponding to the
lation of a surface
performed by the
rdan et al, 2002).
gorithm. For more
sion of its perfor
and above ground
> papers.
ce z(x,y), where
at geographic po-
rary smoothly, the
car combination of
VEI + ly )
V2 + lvyy)]
(2)
Lu, dT, fora
lel can be seen asa
e main philosophy
lered as a surface,
d by buildings and
dd to the high fre-
rain shape would
terms. Hence, the
1e form (2) with a
ound points in the
at is present in the
are estimated from
,m]) of the DEM
obtained from the
dinates into equa-
linear equations
(3)
r(1) SN.NO)
7 (2) SN,N 2)
(m) Sn. nm)
,
th
statistically robust
'stimator theory, in
bove ground DEM
ig a function p of
model predictions
thus is:
js (4)
International Archives of the Photogrammetry, Remote Sensin g and Spatial Information Sciences, Vol XXXV, Part B3. Istanbul 2004
For the function p traditionally Tukey's norm is used. But keep-
ing in mind that ground points systematically yield negative er-
rors co (7), we use the asymmetric adaptation of Tukey's norm
given by:
0 | fe<O
c? € 423 tp
pele) = $ (1-0-1377) ifü«esc
S otherwise
(5)
(cis a scale factor). The weight function w. corresponding to this
norm is:
1 ife
We) = f0ce<e (6)
0 otherwise
Minimization of the object function proceeds iteratively till con-
vergence by the weighted least-squares solution of the (iterated)
system (3) :
OW = (M WPM) M'w? , (7)
where Kk is the iteration step and W^ = diag ( un(e? )) is
the diagonal weight matrix.
The iteration process is started by computing an initial DTM by
means of a least-squares solution ofthe system (3) based on a ran-
dom sample of DEM points from the full DEM, or on the DEM
points contained in the ground level segment as obtained from the
segmentation phase. In the subsequent iteration steps, the value
of c is progressively decreased in order to reject more and more
above ground DEM points as being “outliers”, the weight matrix
w? is updated according to the new value of c and a new model
estimate 8^) is computed. The minimal and the maximal values
of c are user-defined, but they must be chosen in accordance with
the amplitude of the DEM and the minimal height of the above
ground structures in the scene.
4 EXPERIMENTAL RESULTS
The segmentation and the DTM surface fitting algorithm have
been tested on state-of-the-art DEMs obtained through correspon-
dence matching as well as from laser altimetry of several com-
plex urban areas in France with various landscape types ranging
from relatively flat and with modest constant slope to regions with
large variations in terrain slope and altitude. In particular, the fol-
lowing three sorts of data were used in the experiments reported
here.
Synthetic DEM : A synthetic DEM of size 1024 x 1024, which
is generated as follows: First, an analytic “terrain” is computed
as a sum of harmonic functions. Then five “buildings” are added
to the terrain. These are rectangular block shapes with different
heights. Finally, Gaussian noise is added to this analytic scene.
DEM computed by stereo correlation: This DEM is com-
puted with the algorithm of (Cord er al., 2001) on a stereo pair
of scanned aerial images from the Hoengg dataset (Dataset Ho-
engg, 2001). The original images and the DEM are of dimensions
1664 x 1512 and with a ground resolution of about 10 cm per
pixel. The terrain slope on this area is approximately 10 1n.
Airborne laser DEM: The third dataset is an airborne laser
DEM, kindly provided to us by the Institut Géographique Na-
tional (France). The complete DEM covers a 2 km x 2 km area
on the center of Amiens (France), where each pixel has a ground
resolution of 20 em. In the experiments reported here, we used
two sub-areas presenting various terrain slopes: The first one has
dimensions 1212 x 1640 and the second one covers an area of
size 2048 x 2048 in the city center.
First, the segmentation algorithm was applied to the data in order
to test whether large portions of the ground level could be ex-
tracted in a (semi-)automatic manner. It turned out that, regard-
less of the terrain, the same parameter settings (i.e. the choices
for the radii of the spheres — systematically denoted by r in sec-
tion 2, — the scale factor p for scaling the z-coordinates, and
the threshold n for isolated point detection) could be used for
all DEMs with similar dimensions. Due to page restrictions, we
are only able to present here the results obtained from one of the
above mentioned dataset. Readers who are interested in a detailed
discussion of the results of segmentation and DTM estimation on
the other datasets as well are referred to (Van de Woestyne er al.,
2004). It is important to note, however, that the results described
by means of the particular example below transfer to the other
datasets as well.
Figure 3 (a) shows one of the images of the first sub-area of the
Amiens region (France). The corresponding part of the DEM,
which was obtained through airborne laser scanning, is depicted
in Figure 3(b). The coloring in Figure 3 (b) represents the al-
titude of the corresponding scene point (with red indicating the
highest value and blue the lowest). Figure 3(c) illustrates the
effect of smoothing and of isolated point removal on the DEM
in Figure 3(b). The scale factor p applied to the z-coordinates
was set to 10 and the radii of the spheres used for both smooth-
ing and isolated point removal was also set to 10. The threshold
n for an isolated point was set to 4. Observe that mostly DEM
points corresponding to building facades and vegetation are re-
moved. Figure 3 (d) shows the segmentation of the DEM, which
was automatically obtained by the algorithm described in sec-
tion 2.3. Here 12 was used as value for the radius in the definition
of r-connectivity. Remark that a larger value for the radius of the
spheres is used in the segmentation phase than for the preprocess-
ing steps. The reason is that the radius r used in the preprocessing
stage expresses some sort of error tolerance that is applied to the
estimated altitude of the DEM points, whereas the radius r in the
segmentation stage, on the contrary, represents the minimal sep-
aration distance between different surface patches in the scene.
As mentioned at the beginning of section 2, the ground level in
dense urban areas is expected to be largely made up of the road
network; and thus is likely to correspond to one of the largest
DEM segments, to have low altitude when compared to the other
segments, and, most importantly, to extend over the whole urban
area represented in the DEM. Selecting the segment containing
the largest number of DEM points in Figure 3 (d) results in the
area depicted in Figure 3 (e). When a relative altitude coloring is
applied to the segment (i.e. the coloring does not indicate abso-
lute height values, but the colors red and blue respectively corre-
spond to the highest and the lowest point in the segment itself),
then possible errors immediately catch the eye. Indeed, a closer
look to Figure 3(e) shows a red colored spot approximately in
the middle of the figure and near the lower edge, whereas all the
other DEM points in the segment obtain a blueish color. This
indicates that the altitude of the spot seriously deviates from the
average altitude of the rest of the segment. So, probably this is
an error. In the introduction it is mentioned that the segmentation
algorithm is built into a user-friendly, multi-platform software en-
vironment, called ReconLab. ReconLab was initially designed to
