International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol XXXV, Part B3. Istanbul 2004 
yet, and in which all points that are contained within a sphere of 
radius r centered at a segment point are assigned to that segment 
as well. But, the large number of data points in a dense DEM, on 
the one hand side, and the real time requirements we put forward 
for the segmentation algorithm, on the other hand, prompt us for 
another approach. Therefore, segmentation will be performed by 
an indexation procedure, which has the advantage that all DEM 
points have to be visited only once. The procedure comprises the 
following steps. 
1. Start by indexing all the DEM points. If the DEM is con- 
structed over a rectangular regular grid, then the (x, y)-coor- 
dinates of the DEM points can serve as the index. Moreover, 
if grid points would be missing in the DEM (e.g. due to oc- 
clusions in the imagery or because they were omitted during 
isolated point removal), then they may be given an (arbitrary 
negative) z-value which is a few multiples of r smaller than 
the minimal z-value of the DEM points. In this way, it will 
be possible to easily remove them later, while now main- 
taining the regularity of the grid. It is important to note, 
however, that, although a regular grid structure facilitates 
the indexing, it is not required for this algorithm. 
2. Create a segmentation table. This is a dynamic table of lists 
of indices of DEM points that are r-connectable. Initially, 
the segmentation table is empty and will be updated in the 
subsequent steps. 
3. Systematically scan the DEM and perform the following op- 
erations for every DEM point P. 
(a) If the point at hand P does not belong to a previously 
constructed segment, then create a new segment num- 
ber (index) in the segment table and assign this point 
to the new segment. 
Find all DEM points Q that are contained in a sphere 
with radius r centered at P. 
= 
(c) If some of these points Q already belong to a previ- 
ously created segment, then all these segments, to- 
gether with all the other points in the sphere around 
P, are merged by transferring the indexes of the in- 
volved DEM points to the entry corresponding to P 
in the segmentation table (cf. Figure 2). This is be- 
cause all the DEM points in the segment of Q are 
r-connectable to P, and hence, by definition, belong 
to the same segment as P. 
4. Update the segmentation table by removing segment indexes 
that contain no references to DEM points. This can easily 
be performed, because for each segment in the table, the 
number of DEM points in the reference list is stored as a 
separate entry. 
  
  
Segmentation table -. Segmentation table 
  
  
  
  
  
  
  
  
Index | # | References Index | # References 
Q TASSO S qu (pog ih 
Lg i rT Oe ey 
2 |3|0,5,6 3 [310,56 
  
  
  
  
  
  
  
  
  
  
  
  
Figure 2: The segment with index O in the segment table is merged 
with the segment indexed 1 by moving the indexes of the DEM 
points in segment 0 to the reference list of segment 1. 
Observe that, as the number of points in every segment is ex- 
plicitly available, the segments can be sorted accordingly and 
the largest regions (ground level) can easily be identified and se. 
lected. Results obtained by the segmentation algorithm on state- 
of-the-art DEMs are presented and discussed in section 4. 
3 DTM ESTIMATION 
The DEM points that belong to the segment corresponding to the 
ground level serve as input data for the estimation of a surface 
model for the terrain. The DTM estimation is performed by the 
algorithm presented in (Belli ef al., 2001, Jordan et al., 2002), 
This section provides a brief overview of this algorithm. For more 
details on the estimation method and a discussion of its perfor- 
mance on a full DEM (i.e. including buildings and above ground 
areas), the interested reader is referred to those papers. 
A DTM is modeled as a parameterized surface z(x,y), where 
z(x,y) denotes the altitude of the scene point at geographic po- 
sition (z, y). As the terrain is assumed to vary smoothly, the 
function z(x, y) can be decomposed into a linear combination of 
2D harmonic functions: 
z(x,y) = aoo 
N 
num. 
k.lzzO;k-Figz0 
N 
bari, 
K,1z:0 ; k-4-0:50 
ar1 cos (2x (kv2x + lvyy)] 
bu sin (2x (kvax + lvyy)] 
(2) 
with fundamental frequencies #, = 1/T, and v, = 1/T, fora 
DTM of size T. x T,,. The order N of the model can be seen asa 
constraint on the terrain variability. In fact, the main philosophy 
behind the method is that, if a DEM is considered as a surface, 
then the variations in altitude (z-values) caused by buildings and 
other above ground structures would mainly add to the high fre- 
quencies of the decomposition, whereas the terrain shape would 
mainly be represented by the low frequency terms. Hence, the 
idea is to model the DTM by a function of the form (2) with a 
low value for N, and to consider the above ground points in the 
DEM as being "structural noise" (outliers) that is present in the 
data. 
The 2(N 4- 1)? — 1 parameters in equation (2) are estimated from 
the coordinates (z;, yi, z;) (à € (1,2,..., m]) of the DEM 
points contained in the ground level segment obtained from the 
segmentation phase. Substituting these coordinates into equa- 
tion (2) leads to an overdetermined system of linear equations 
z—-MO (3) 
ith z — z zm) 
wit z «(201.2204 Am), 
1 Cp,1(1) 89,1 (1) CN,NU) Sx nN) 
Co,1(2) So.1(2) CN.N (2) SN,N@) 
M = 
1 €. 1m) So.1 On, ) CN,N (M) Sn. nm) 
where Cy, 1(i) — cos [2m (kve xi - Uv, yi) ]. 
Ska (i) = sin [27 (k vs x; + lvy y:) | and with 
© = (ao,0; 40,1; bo,1;-- - ; AN,N, bN,n)". A statistically robust 
solution to system (3) is obtained by using M-estimator theory, in 
which the influence of outliers — in this case, above ground DEM 
points — is iteratively reduced by minimizing a function p of 
the errors ce (i) along the z-axis between the model predictions 
M®(i) and the data z; . The optimal solution thus is: 
m 
©, = arg min > pleali)y) . (4) 
$m 
  
International Ai 
EE t 
For the functior 
ing in mind tha 
rors co (i), we 
given by: 
pele) = 
(cis a scale fact 
norm is: 
we(e) = 
Minimization of 
vergence by the 
system (3) : 
^ 
eU 
where k is the if 
the diagonal wei 
The iteration pr 
means of a least- 
dom sample of | 
points contained 
segmentation ph 
of c is progressi 
above ground D 
w^ is updated 
estimate (^ is 
of c are user-def 
the amplitude of 
ground structure 
The segmentatic 
been tested on st: 
dence matching 
plex urban areas 
from relatively fl 
large variations i 
lowing three sor 
here. 
Synthetic DEM 
is generated as f 
as a sum of harm 
to the terrain. TI 
heights. Finally, 
DEM computec 
puted with the a 
of scanned aerial 
engg, 2001). The 
1664 x 1512 an 
pixel. The terrair